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PREFACE 


The papers included herein were presented at the NASA Tire Modeling Workshop 
held at Langley Research Center, September 7-9, 1982. The workshop was organized 
into six sessions dealing with finite element developments, applications to tire 
dynamic problems, solution techniques for tire contact problems, experimental 
data, tire thermal studies, and current design practices. Discussion periods at 
the end of each session allowed the participants to describe their own work or 
problem areas and to make recommendations for future directions of tire modeling 
research. 

The objective of the workshop was to provide a forum for the interchange of 
information among tire analysts and to establish goals and objectives for future 
research programs in this area. Hence the various conclusions and recommendations 
suggested by the workshop participants are expected to provide guidance for tire 
modeling studies within NASA and, to some extent, the tire industry for the next 
several years. 

Special recognition is due Ms. Connie 0. Featherston for her planning, 
coordination, and execution of the various administrative tasks associated with 
the workshop. 

The use of trade names or names of manufacturers in this report does not 
constitute an official endorsement of such products or manufacturers, either 
expressed or implied, by the National Aeronautics and Space Administration. 


John A. Tanner 
Langley Research Center 
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FINITE ELEMENT MODELING AND ANALYSIS OF TIRES 


Ahmed K. Noor 

Joint Institute for Advancement of Flight Sciences 
George Washington University 

and 

C. M. Andersen 
College of William and Mary 


Abstract 

Although the problem of tire modeling and analysis has been a subject of 
continuing concern for the tire industry, to date no simple and general tire model 
exists for predicting the response of the tire under various loading conditions. 

Much of the recent progress in finite element technology has not been exploited for 
tire modeling and analysis. The present paper focuses on this issue. Specifically, 
the paper reviews some of the recent advances in finite element technology which 
have high potential for application to tire modeling problems. It also identifies 
the analysis and modeling needs for tires. 

The topics covered include: 1) reduction methods for large-scale nonlinear 

analysis, with particular emphasis on treatment of combined loads, displacement- 
dependent and nonconservative loadings; 2) development of simple and efficient 
mixed finite element models for shell analysis, identification of equivalent mixed 
and purely displacement models, and determination of the advantages of using 
mixed models; and 3) effective computational models for large-rotation nonlinear 
problems, based on a total Lagrangian description of the deformation. 


INTRODUCTION 


The problem of tire modeling and analysis has long been an area of major con- 
cern to the tire and aircraft industries. A hierarchy of models varying in the 
degree of sophistication has been proposed. Some of these models are listed in 
Fig. 1 and are sketched in Fig. 2. For a detailed description of the models see 
Ref. 1. The models are grouped into six groups as follows: 

The first group consists of the early tire models which are characterized by 
their simplicity. Among these models are the ■6tAA.ng, b&am, and fhLng on ztcu,td.c. {oH. 
V^COeXaS-tZc) foundations. These models were used by Clark and co-workers (Ref. 2). 
Their major drawbacks are: 1) they require extensive experiments to evaluate 
the equivalent properties, and 2) their accuracy and range of validity are not 
known in advance. 

The second group consists of the c.Oh.d-noJMOhh modol^ , which are sometimes 
referred to as netting anaJLysiS , wherein the inflation pressure is assumed to be 
carried exclusively by the cords (see Ref. 3) . These models have the drawback of 
neglecting both the bending in the tire and the stiffening effect of the rubber. 

The third group of models are the meJvbAane. mod^is, which are based on the use 
of a linear or nonlinear momentless theory of shells (Refs. 4, 5 and 6). Their 
major drawback is that they cannot handle discontinuities in loading, geometry or 
material properties. 

The fourth group is the two-dime^niionoZ axLsymm^t^c, modoLs (Ref. 7), which 
are limited to axlsymmetric loadings. 

The fifth group is the thfLZQ,-cUjne,n6ionat aontinuum modo.Zs. Two approaches 
have been proposed for the analysis of these models. The first approach is based 
on using seml-analytlc techniques to reduce the dimensionality of the problem 
(e.g., Fourier expansions in the circumferential direction). The second approach 
is based on using three-dimensional isoparametric solid elements. 

The sixth group of models Includes a variety of tu)0- dimensional, thin and thick 
shell models (see, for example. Refs. 8 and 9). Thin shell models neglect trans- 
verse shear deformation, and their use for modeling tires is therefore 
questionable. Anisotropy results in increasing the size of the analysis model, 
and consequently many investigators neglect its effects by using an orthotropic 
model. 

The present paper focuses on the use of two-dimensional thick shell models. 


2 



IIFFEIEIIT HIE MOBELS 


MODEL 

DRAWBACKS 

• SIMPLE MODELS 

STRING) ELASTIC 

BEAM ON OR 

RING 1 VISCO- 

ELASTIC 

FOUNDATION 

• REQUIRE EXTENSIVE EXPERIMENTS TO 
EVALUATE EQUIVALENT PROPERTIES 

• ACCURACY AND RANGE OF VALIDITY 
NOT KNOWN IN ADVANCE 

• CORD-NETWORK MODELS 

• NEGLECT BENDING IN TIRE AND 
STIFFENING EFFECTS OF RUBBER 

• MEMBRANE MODELS (BASED ON 
MOHENTLESS THEORY OF SHELLS) 

• CANNOT HANDLE DISCONTINUITIES (OR 
SHARP CHANGES) IN LOADING, GEOMETRY 
(CURVATURE) OR MATERIAL PROPERTIES 

• TWO-DIMENSIONAL AXISYMMETRIC MODELS 

• LIMITED TO AXISYMMETRIC LOADING 

• THREE-DIMENSIONAL MODELS 

• SEMI-ANALYTIC SOLUTIONS 

• THREE-DIMENSIONAL SOLID ELEMENTS 

• COMPUTATIONALLY EXPENSIVE 

• LAMINATED ANISOTROPIC SHELL MODELS 
(THIN AND THICK) 

• THIN SHELL MODELS NOT ADEQUATE 

• EFFECT OF ANISOTROPY CAN BE 
SIGNIFICANT 


Figure 1 



Figure 2 
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TIRE CONFIGURATION AND COMPONENTS 


T3^ical configurations and components of modern tires are shown in Fig. 3 (see 
Ref. 10). Commercially successful tires are now built as a series of layers of 
flexible high-modulus cords encased in a low-modulus rubber or rubber-like material. 
Hence, a laminated (or layered) model is needed. 



I 


TIRE LOADS 

The three types of loads applied to the tire and their major characteristics 
are listed in Fig. 4. The three load tjrpes are: 

1) In^i/ltion p-'teAAUA.e, which is axisymmetric but is displacement dependent 

2) MzakayiLcaJi Zoadi which include centrifugal force. Impact loading, contact 
forces, and frictional forces; except for the centrifugal force, which is axisym- 
metric (and displacement dependent) , all the other loads are symmetric 

3) TheAmaZ ZoacU, which arise due to various manufacturing and operating 
conditions, such as unequal expansion and contraction of rubber and cord, 
hysteretic heating, sliding of the tread on a rough surface, and cord shrinking 
after molding 


LOADS 

CHARACTERISTICS 

• INFLATION PRESSURE 

• AXISYMMETRIC BUT DISPLACEMENT- 
DEPENDENT 

• MECHANICAL LOADS 

• CENTRIFUGAL FORCE 

• IMPACT LOADING 

• CONTACT FORCES 

• FRICTIONAL FORCES 

• AXISYMMETRIC BUT DISPLACEMENT- 
DEPENDENT 

• ASYMMETRIC 

• 

• 

1 

1 

• THERMAL LOADS 

• UNEQUAL EXPANSION AND 
CONTRACTION OF RUBBER 
AND CHORD 

• HYSTERETIC HEATING 

• CHORD SHRINKING AFTER 
MOLDING 

• ASYMMETRIC 


Figure 4 
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CHARACTERISTICS OF EFFECTIVE SHELL ELEMENTS 
FOR ANALYZING TIRES 

The characteristics of an effective shell finite element model for analyzing 
tires are listed in Fig. 5. The shell element is developed using either a con- 
sistent two-dimensional shell theory or a three-dimensional continuum theory with 
proper interpolation functions in the thickness direction. The elements obtained by 
using the latter approach are referred to as de.gzn2Aatz 6h&Zt eZe.me.yiti, If a two- 
dimensional shell theory is used, the elements need to be deep and curved and must 
account for each of the following effects: 

1) Laminated construction and anisotropic material behavior 

2) Variation in geometry (e.g., curvature and thickness) as well as of other 

lamination parameters 

3) Transverse shear deformation 

4) Large rotations 

5) Pressure stiffness (for displacement-dependent loadings such as inflation 

pressure) 

6) Thermovlscoelastic material response 


• BASED ON EITHER 

• CONSISTENT TWO-DIMENSIONAL SHELL THEORY, OR 

• THREE-DIMENSIONAL CONTINUUM THEORY WITH PROPER INTERPOLATION FUNCTIONS 
IN THE THICKNESS DIRECTION (DEGENERATE SHELL ELEMENTS) 

• DEEP, CURVED ELEMENTS 

• INCLUDE EFFECTS OF: 

• LAMINATED CONSTRUCTION AND ANISOTROPIC MATERIAL BEHAVIOR 

• VARIATION IN GEOMETRY (E.G., CURVATURE AND THICKNESS), LAMINATION PARAMETERS 

• TRANSVERSE SHEAR DEFORMATION 

• LARGE ROTATIONS 

• PRESSURE STIFFNESS (FOR DISPLACEMENT-DEPENDENT LOADING; E.G., INFLATION 
PRESSURE) 

• THERMOVISCOELASTIC MATERIAL RESPONSE 


Figure 5 
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OBJECTIVES AND SCOPE 

The objectives of this paper are listed In Fig. 6. They are: 

1) To review some recent developments In finite element technology which are 

applicable to the analysis and modeling of tires 

2) To Identify some of the analysis and modeling needs for tires 

The paper Is divided into four parts. The first part deals with new develop- 
ments In reduction methods for nonlinear problems. These Include computational 
procedures for handling combined, displacement-dependent, and nonconservative 
loads. The second part of the paper deals with mixed finite element models for 
tires in which the fundamental unknowns consist of both force and displacement 
parameters. The equivalence of some of these models with some of the purely 
displacement models is discussed. ' 

The third part of the paper deals with large-rotation nonlinear problems. Two 
formulations are presented; namely, a mixed formulation and a penalty formulation. 
Both formulations are based on the total Lagrangian description of the deformation. 
The fourth and last part of the paper deals with analysis and modeling needs for 
tires. 


OBJECTIVES 

t REVIEW SOME RECENT DEVELOPMENTS IN FINITE ELEMENT TECHNOLOGY WHICH ARE 
APPLICABLE TO ANALYSIS AND MODELING OF TIRES 

• IDENTIFY ANALYSIS AND MODELING NEEDS FOR TIRES 
SCOPE 

• REDUCTION METHODS FOR NONLINEAR PROBLEMS 

• COMBINED LOADING PROBLEMS 

• DISPLACEMENT-DEPENDENT AND NONCONSERVATIVE LOADING PROBLEMS 

• DYNAMIC PROBLEMS 

• MIXED FINITE ELEMENT MODELS 

• EFFICIENT AND ACCURATE MIXED MODELS 

t EQUIVALENT CLASSES OF MIXED MODELS AND REDUCED/SELECTIVE INTEGRATION 
DISPLACEMENT MODELS 

• MERITS OF MIXED MODELS OVER EQUIVALENT DISPLACEMENT MODELS 

• LARGE-ROTATION NONLINEAR PROBLEMS 

• MIXED FORMULATION 

• PENALTY FORMULATION 


Figure 6 
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REDUCTION METHODS FOR NONLINEAR TIRE PROBLEMS 


The first topic considered in this paper is reduction methods for nonlinear 
analysis. The basic features of reduction methods are outlined in Fig. 7. They 
are techniques for reducing the number of degrees of freedom through the trans- 
formation shown in the figure. The vector {X} represents the original displacement 
degrees of freedom. The vector refers to amplitudes of displacement modes and 
[r] is a transformation matrix whose columns represent a priori chosen global 
displacement modes. 


As is to be expected, the effectiveness of reduction methods depends to a great 
extent on the proper selection of the displacement modes. In a number of studies 
it was shown that an effective choice of the displacement modes includes the various- 
order derivatives of the displacement vector with respect to the load parameter (see 
Refs. 11 and 12). These vectors are gmeJmt.e.d by ii6^ng the. eZe.me.nt modeZ 


2 

8X 9 X 

the. ttAe., The recursion formulas for evaluating the derivatives {r— }, { — ;r}, ... 

are obtained by successive differentiation of the original finite element equations. 
The left-hand sides of the recursion formulas are the same (see Ref. 12). Therefore 


onZy one. matnZx {^actonZzatton Jj> n.equuA.ed ^oh. the genenxxtZon aZZ the gtoboZ 
appflOXAJnatZon vector. Several numerical experiments have demonstrated the effect- 
iveness of this choice (see Refs. 12 and 13). 


DEFINITION: ARE TECHNIQUES FOR REDUCING THE NUMBER OF D.O.F. 
THROUGH THE TRANSFORMATION 

|X| = ORIGINAL D.O.F. IN THE FINITE ELEMENT 
MODEL 

[r] = MATRIX OF GLOBAL DISPLACEMENT MODES 

= REDUCED D.O.F. -AMPLITUDES OF DISPLACEMENT. 
MODES 


JUSTIFICATION: FOR MANY TIRE PROBLEMS THE LARGE NUMBER OF 

D.O.F. IXt IS DICTATED BY THE COMPLEX TOPOLOGY 
OF THE TIRE (DISCONTINUITIES IN GEO/VETRY, 
LAMINATION, ETC.) RATHER THAN BY EXPECTED 
COMPLEXITY OF BEHAVIOR 


SELECTION OF GLOBAL DISPLACEMENT MODES : 


= 

■jdXj (d2X| 

1 

d^XI 

■ 


Lp'l 

1 

r 


= 

LOAD PARAMETER 



• COLUMNS OFCn GENERATED BY USING THE ORIGINAL FINITE ELEMENT 
MODEL OF THE TIRE 

• THEIR GENERATION REQUIRES ONLY ONE 'LARGE MATRIX' FACTORIZATION 

• NUMERICAL EXPERIMENTS HAVE DEMONSTRATED THEI R EFFECTIVENESS 


Figure 7 
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BASIC EQUATIONS USED IN REDUCTION jMETHODS 
FOR NONLINEAR TIRE PROBLEMS 


The basic equations used in the reduction methods for geometrically nonlinear 
tire problems are given in Fig. 8. It is worth noting that the original displace- 
ment unknowns {X} can be on the order of thousands whereas the reduced unknowns {i(j} 
are typically 20 or less. This is true regardless of the complexity of the structure 
and/or loading. The details of the computational procedure for tracing the load- 
deflection paths in geometrically nonlinear static analysis are given in Refs. 11 
and 12. 



ACTUAL (LARGE) PROBLEM 

REDUCED (SMALL) PROBLEM 

FUNDAMENTAL 

UNKNOWNS 



1X(= INDIVIDUAL - - AMPLITUDES 

DISPLACEMENTS 
• THOUSANDS OF UNKNOWNS 

Uh UlbFLALtmtNl 

MODES 

• TWENTY OR LESS 

GOVERNING 

EQUATIONS 

[KlIXt + IG (X)}- piPl = 0 
~ 1000 EQUATIONS 

[K]{(AI + IG ((^)f-plPl = 0 
~ 20 EQUATIONS 

HOW TO TRACE LOAD- 
DEFLECTION PATH 

• REPEATED SOLUTION OF 
LARGE SYSTEMS OF 
SIMULTANEOUS NONLINEAR 
ALGEBRAIC EQUATIONS 

• GENERATION OFCH 

• MARCHING WITH SMALL 
SYSTEM OF EQUATIONS 

• ERROR SENSING AND 
CONTROL (UPDATING [T] 
WHENEVER NEEDED) 


[K] = cr]T[K][r] . \G(iif)\ =Ln'^\Gm\ . iPi = cr]TlPl 


Figure 8 
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APPLICATION OF REDUCTION METHODS TO GEOMETRICALLY NONLINEAR 
ANALYSIS OF A TIRE SUBJECTED TO UNIFORM INTERNAL PRESSURE 

As a simple application of reduction methods to the geometrically nonlinear 
analysis of tires, consider the laminated anisotropic elliptic toroidal shell 
shown in Fig. 9. Due to axial symmetry, only one meridian was modeled using 
four-noded elements with cubic Lagrangian interpolation functions for all the 
displacement and rotation degrees of freedom. The high accuracy of the total 
strain energy obtained by using six basis vectors is demonstrated in Fig. 9. 



Gjj = 270 psl h = 0,1)2 In. 

>'l_j =0.4 e = 0.5 


CORD ORIENTATION + 45/- 45/ + 45/ - 45 ... . 

NL = 10 

Figure 9 
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TREATMENT OF COMBINED LOADS 


The basic equations and the computational procedure used in applying reduction 
methods to the analysis of tires subjected to combined loads are highlighted in 
Figs. 10 and 11. For simplicity, only two independent loads are considered. 

First, the original finite element equations are given. The external loading 
is normalized with respect to two independent load parameters p^ and P 2 ■ The basis 
reduction is done as before, via the transformation shown in Fig. 10. Then the 
Raylelgh-Ritz technique is used to approximate the original set of finite element 
equations by a reduced system of equations in the new unknown parameters {if)}. The 
number of these equations is considerably less than that of the original equations. 

As previously noted, the crux of reduction methods is the proper selection of 
the transformation matrix [r]. In the case of combined loading, the columns of the 
matrix [p] are selected to be the various-order derivatives of the displacement 
vector {x} with respect to the two independent parameters p^ and 

To trace the different nonlinear paths, corresponding to different combinations 
of the independent load parameters, the basis vectors are evaluated for the unloaded 
structure (p-]^ = P 2 = 0) , and the corresponding reduced equations are generated. The 
different nonlinear paths of the tire are obtained by fixing one of the load para- 
meters, varying the other, and repeating the process with different values of the 
first load parameter. This is all done using the -SOJne, iQJt fLZducQ-d ZqiMltdoni , 

The total cost of the analysis, to a first approximation, is little more than the 
cost of one linear solution of the original, full system of finite element equations. 
The procedure is described in detail in Ref. 14. 

As a by-product of this technique, a considerable reduction can be made in the 
size of the analysis model used in studying the nonlinear response of tires sub- 
jected to asymmetric loading. This can be accomplished by decomposing the loading 
into S3mimetrlc and antls 3 rmmetric components and treating each as an independent 
loading. 


REDUCTION METHODS FOR NONIINERR PROOIEMS 
TREATMENT OF COMOINEO LOADS 


GOVERNW6 FINITE ELEMENT EOUATIONS 


[KHx} + {G(X)} - = 0 

[K] = LINEAR GLOBAL STIFFNESS MATRIX 
{X} = VECTOR OF NODAL DISPLACEMENTS 
{G(X)} = VECTOR OF NONLINEAR TERMS 

= NORMALIZED LOAD VECTORS 
Pj^, P 2 = independent LOAD PARAMETERS 

BASIS REDUCTION 




p 


2 


{X} 


n, 1 


[F] 


iT,r^'''Jr,l 


r << n 


REDUCED SYSTEM OF EQUATIONS 

+ (Gd^)} - p^{p'^>} - p^TP^^’} = 0 


Figure 10 


TREATMENT OF COMBINED LOADS 

SELECTION OF BASIS VECTORS 


{X} = [rirt) 


[ri = 


■3Pi 




3PI 




) { 


3p 


2 / • • • 

3Po 


COMPUTATIONAL PROCEDURE 


• EVALUATE BASIS VECTORS AT Pj = P 2 = 0 (UNLOADED TIRE) AND GENERATE 
REDUCED EQUATIONS 

• TRACE DIFFERENT EQUILIBRIUM PATH BY FIXING ONE OF THE LOAD 
PARAMETERS AND VARYING THE OTHER (USING THE SAME SET_OF 
REDUCED EQUATIONS) 




P 2 


X 




NOTE ; THIS APPROACH CAN BE USED TO REDUCE THE SIZE 

OF ANALYSIS MODELS FOR THE CASE OF UNSYMMETRIC LOADINGS. 


1 2 


Figure 11 
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TREATMENT OF DISPLACEMENT-DEPENDENT 
AND NONCONSERVATIVE LOADING 

The basic equations used In applying reduction methods to the nonlinear analysis 
of tires subjected to displacement-dependent loading are given In Fig. 12. 

First, the governing finite element equations for the total Lagranglan formu- 
lation are shown. The only new term In these equations Is the pressure stiffness 
matrix which represents the follower-load effect, and Is unsymmetrlc for noncon- 
servative loadings. The basis reduction Is done and the reduced equations are 
obtained In the manner outlined previously. The following two Important facts are 
to be noted: 

1. The basis vectors are evaluated for the unloaded structure. Hence, the 

pressure stiffness matrix does not enter Into the left-hand side and onty the. LLneoJi 

6ymmetJUa gtobaZ 6tt^^ne^-6 matJu.x med& to be deaompo6&d. 

2. Since the reduced equations are small In number (on the order of ten or 

less) no S5rmmetrization is needed in the case of nonconservative loadings. 


GOVERNING FINITE ELEMENT EQUATIONS 

FOR A TOTAL LAGRANGIAN FORMULATION 

[[K] - p[K^P^]] {X} + {G(X)> - p{P} = 0 

= PRESSURE STIFFNESS MATRIX (UNSYMMETRIC FOR 
NONCONSERVATIVE LOADING) 

BASIS REDUCTION 

' {x} = [r] 

REDUCED SYSTEM O F EQUATIONS 

|[K] - p[K^Ph]{i(;} + {G(ip)} - p{P} = 0 

NOTES : 

• BASIS VECTORS ARE EVALUATED AT p = 0. THEREFORE, ONLY THE SYMMETRIC [k] 
MATRIX NEEDS TO BE DECOMPOSED. 

• REDUCED EQUATIONS ARE SMALL IN NUMBER (~10). THEREFORE, NO SYMMETRIZATION 
IS NEEDED, 

Figure 12 
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APPLICATION OF REDUCTION METHODS TO THE BIFURCATION BUCKLING 
ANALYSIS OF A RING SUBJECTED TO HYDROSTATIC PRESSURE 


As a simple application of reduction methods to structures subjected to dis- 
placement-dependent loadings, consider the circular ring subjected to hydrostatic 
pressure shown in Fig. 13. 

Doubly-symmetrlc buckling modes are considered; hence, only one quadrant of the 
ring was analyzed using higher-order shear-flexible elements with a total of 59 non- 
zero degrees of freedom. The lowest three buckling loads obtained using three, 
four and five vectors are listed in Fig. 13. The lowest buckling load obtained by 
using four vectors agrees, to five significant digits, with that obtained using the 
full system of equations. With five vectors, the error in the third buckling load 
is less than 3%. 


NUMBER OF 
BASIS 

EIGENVALUES p 

P 

El 

VECTORS 

Pi 

h 

> 1 . 

p? 

3 

2. 9998 

15.1888 


4 

2.9997 

14.9964 

36.5126 

5 

2.9997 

14.9957 

35.9668 

FULL SYSTEM 
(59 D.O.F.) 

2.9997 

14.9947 

34.9751 
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REDUCTION METHODS FOR NONLINEAR DYNAMIC PROBLEMS 


The application of reduction methods to transient tire problem is highlighted 
in Fig. 14. First, the governing semi-discrete finite element equations are given 
for the case of no damping. Then the key elements for an effective reduction 
method are listed (see Ref. 13). They include: 

1) The proper selection of basis vectors (the columns of the matrix [p]) 

2) Characterization of nonlinear dynamic response by means of one or few 

scalars 

3) Sensing and controlling the error in the reduced system of equations 


GOVERNING SEM I- DISCRETE PtNITE ELEMENT EQUATIONS 
fM) {XI = |P}^ - [K] {XI - {G(X)| 


KEY ELEMENTS FOR EFFECTIVE REDUCTION METHOD 


' PROPER SELECTION OF BASIS VECTORS 

CHARACTERIZATION OF NONLINEAR DYNAMIC RESPONSE BY 
MEANS OF ONE OR FEW SCALARS 


• SENSING AND CONTROLLING THE ERROR IN THE REDUCED 
SYSTEM OF EQUATIONS 


Figure 14 
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SELECTION OF BASIS VECTORS FOR THE CASE OF STEP LOADING 


A particular choice of basis vectors which was found to work well for the case 
of step loading is shown in Fig. 15. The vectors consist of a few eigenvectors of 
the linear problem and a few eigenvectors of the steady-state (static) nonlinear 

9G. 

problem. The matrix [„ ■ ] is obtained by using the steady-state (static) nonlinear 

3 

solution. Reduction methods can be used to reduce the computational effort required 
for generating the steady-state nonlinear solution. 


P(t) 


BASIS VECTORS CONSIST OF: 

• FEW EIGENVECTORS OF LINEAR PROBLEM 

[K1 {XI = XIMl {X} TIME t 


• FEW EIGENVECTORS OF STEADY-STATE (STATIC) NONLINEAR PROBLEM 


(KI 


dGj^ 


dX. 

J 


IXI 


XfMl IXl 


Figure 15 
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APPLICATION TO CLAMPED SHALLOW SPHERICAL CAP SUBJECTED 
TO A CONCENTRATED LOAD AT THE APEX 


As a simple application of reduction methods to nonlinear dynamic problems, 
consider the clamped spherical cap subjected to a concentrated load which has a step 
variation in time (see Fig. 16) . The displacement time history obtained using the 
full system of finite element equations, the reduced system with ten linear vibra- 
tion modes, and the reduced system with the proposed set of modes are shown in 
Fig. 16. Tho. bo6^ ueoioAA (eZgzmode^) weAe. not apdatzcl thAoughout the. anaty^tA. 

As can be seen from Fig. 16, the proposed set of basis vectors predicts qualitatively 
the correct response. The phase shift was almost eliminated by increasing the 
number of basis vectors to 14. 



TIME (MICROSECONDS) 

Figure 16 


17 





MIXED FINITE ELEMENT MODELS FOR TIRES 


The second topic considered in this paper is mixed finite element models for 
tires. The basic features of the mixed models are outlined in Fig. 17. The finite 
element models include the effects of both laminated anisotropic construction and 
transverse shear deformation, and allow the geometric and material properties to 
vary within individual elements. The fundamental unknowns consist of the eight 
stress resultants and the five generalized displacements. The stress resultants 
are discontinuous at element Interfaces, and therefore can be eliminated on the 
element level. 


• TIRE MODELED USING LAMINATED ANISOTROPIC. SHEAR-FLEXIBLE. DEEP SHELL 
ELEMENTS WITH VARIABLE GEOMETRIC AND MATERIAL PROPERTIES 

• FUNDAMENTAL UNKNOWNS ARE; 

STRESS RESULTANTS N„g . M„g . Q„ 

GENERALIZED DISPLACEMENTS . w . 

• STRESS RESULTANTS ARE DISCONTINUOUS AT ELEMENT INTERFACES - ELIMINATED 
ON ELEMENT LEVEL 
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MATHEMATICAL FORMULATION FOR THE MIXED MODEL 

The mathematical formulation for the two-field mixed model is based on the use 
of a moderate-rotation nonlinear shell theory in conjunction with the Hellinger- 
Reissner mixed variational principle. The basic features of this formulation are 
outlined in Fig. 18. Different approximation functions are used for each of the 
stress-resultant fields and the generalized displacement field. The governing finite 
element equations for individual elements can be partitioned as shown in Fig. 18. 

The vector {M} is quadratic in {X} and the vector {G} is bilinear in {H} and {X}. 

For mixed models with discontinuous stress resultants at element interfaces, 
the stress resultants can be eliminated on the element level and the governing finite 
element equations reduce to cubic equations in {X} (see Ref. 15). 

APPROXIMATION FUNCTIONS 
STRESS RESULTANTS 


= WUH} . |H|= VECTOR OF STRESS N 

RESULTANT PARAMETERS 

DISPLACEMENTS 


= [yV]|X| . I X}= VECTOR OF NODAL 

DISPUCEMENTS 

G OVERNIN G FINITE ELEMEN T EQ UATIONS FOl 

■ -[F][s]i jwm 

,[S]T [OjJ\xj V|G(H.X)tl IpJ 
WHERE |/M(X)| ANDIGIH, X)| ARE VECTORS OF NONLINEAR TERMS 
(QUADRATIC AND BILINEAR IN{H}AND|X{). 

DISCONTINU OUS STR ES S RESUL TA NTS AT ELEMENT INTERFACES 
Ih] = [F]-l [S]{X} + [Fr^ |/W(X)} 

AND GOVERNING FINITE ELEMENT EQUATIONS REDUCE TO 
[S]T [F] -1 [SJIXI + |G(X)} = IP} 

WHERE {G(X)| = VECTOR OF NONLINEAR TERMS (CUBIC IN )Xl). 


Figure 18 
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REDUCED /SELECTIVE INTEGRATION DISPLACEMENT MODELS 


In recent years a class of displacement models with a performance comparable 
to that of the mixed model has been developed. These are referred to as reduced/ 
selective integration displacement models (see, for example. Refs. 15 and 16). 

The major features of these models are outlined in Fig. 19. The governing finite 
element equations for the individual elements are cubic in {X}. The definitions 
of full, reduced, and selective integration are given in Fig, 19. If an nxn Gauss- 
Legendre formula is used to Integrate the linear stiffness matrix [k] exactly for 
parallelogram elements, then in full integration nxn quadrature points are used. 

In reduced integration (n-1) x (n-1) quadrature points are used, and in selective 
integration nxn quadrature points are used for some terms of [k] and (n-1) x (n-1) 
points for other terms. 


GOVERNING FINITE EliMENT EQUATIONS FOR INDIVIDUAL ELEMENTS 
[KliXl + |G(X)| = |Pl 

{X} = VECTOR OF NODAL DISPLACEMENTS 
|G(X){ = VECTOR OF NONLINEAR TERMS (CUBIC IN {X{) 

FULL (NORMAL). REDUCED AND SELECTIVE INTEGRATION 

• IF nxn GAUSS-LEGENDRE FORMULA IS USED TO INTEGRATE [K] EXACTLY 
FOR RECTANGULAR (OR PARALLELOGRAM) ELEMENTS 

• FULL (NORMAL) INTEGRATION USES nxn QUADRATURE POINTS 

• REDUCED INTEGRAWN USES (n-1) x (n-1) QUADRATURE POINTS 

• SELECTIVE INTEGRAWN USES nxn POINTS FOR SOME TERMS AND 
(n-1) X (n-1) POINTS FOR OTHER TERMS OF [K] 


Figure 19 



EQUIVALENT FINITE ELEMENT MODELS 


The equivalence between finite element models is defined in Fig. 20. Finite 
elements are equivalent if their individual governing equations, when expressed in 
terms of a common set of nodal variables and/or parameters, are identical (see 
Ref. 16). It is important to note that the other parameters not contained in the 
common set are local to the individual elements. For nearly equivalent models, 
the finite element equations are almost identical. 


DEFINITION OF EQUIVALENCE 

• WO FINIII ELEMENT MODELS ARE EQUIVALENT IF THEIR GOVERNING 
FINITE ELEMENT EQUATIONS, WHEN EXPRESSED IN TERMS OF A COMMON 
SET OF NODAL VARIABLES AND/OR PARAMETERS, ARE IDENTICAL . 

THE OTHER NODAL VARIABLES AND/OR PARAMETERS NOT CONTAINED 
IN THE COMMON SET MUST BE LOCAL TO THE INDIVIDUAL ELEMENTS 
(I.E., DO NOT AFFECT THE ASSEMBLY PROCESS). 

• NEARLY EQUIVALENT MODELS ARE ONES FOR WHICH THE GOVERNING 
FINITE ELEMENT EQUATIONS ARE ALMOST IDENTICAL. 


Figure 20 
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EQUIVALENT MIXED AND DISPLACEMENT MODELS 


The governing finite element equations of both mixed and displacement models 
are shown in Fig. 21 and the mathematical requirements for the equivalence of the 
two models are listed. The table given in Fig. 21 lists examples of equivalent 
quadrilateral (in planform) mixed and displacement models. The following can be 
noted: 

1) Equivalent mixed and displacement models have the same number of displace 
ment nodes, use the same approximation functions for the generalized displacements 
and use the same number of numerical quadrature points. 

2) The sjmibol (F) refers to full integration and (R) refers to reduced 
integration. 

3) If the geometric and material characteristics within the individual 
elements are constants, the number of quadrature points listed in the table 
generates exact integrals for the mixed models and only approximate Integrals 
for the displacement models. 


EQ UIVALENCE 


MIXED MODEL WITH DISCONTINUOUS 
STRESS RESULTANTS 

DISPLACEMENT MODEL 

GOVERNING FINITE ELEMENT EQUATIONS 


[S]^[F]'^ [S]{X}+{6(Xt! = IP} 

[K]|X} + |G(X)|= jPl 

FROM WHICH 


[S]^ [F]'^CS] = 

= [K] 

ie(X)i = 

= |G(X)1 


NEAR EQUIVALENCE = IS REPLACED BY - 
EXAMPLES OF EQUIVALENT QUADRILATERAL ELEMENTS 


NUMBER OF 
DISPLACEMENT 
NODES 

QUADRATURE 

POINTS 

NUMBER OF STRESS 
RESULTANT PARAMETERS 


4 

2 X 2 (F) 
1 (R) 

4 

1 

2:^ 

9 

3x 3 (F) 
2 X 2 IR) 

9 

4 


16 

4 X 4 (F) 
3 X 3 (R) 

16 

9 

0 


Figure 21 
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NONLINEAR RESPONSE OF CIRCULAR TOROIDAL SHELL 
SUBJECTED TO UNIFORM EXTERNAL PRESSURE 


To assess the accuracy of the mixed models with discontinuous stress resultants 
at interelement boundaries, the nonlinear response of the circular toroidal shell 
shown in Fig. 22 is analyzed using these models. The solutions obtained using six 
and eight finite elements with nine displacement nodes and four stress nodes are 
compared with the converged solution in Fig. 22. 



E = 1 X 10^ psl 
p = 100 psl 
a = 15. In. 
b = 10. In. 



X 


Figure 22 
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CLAMPED CYLINDRICAL SHELLS SUBJECTED 
TO UNIFORM PRESSURE LOADING 

To assess the accuracy of the different displacement and mixed models, the 
large-deflection nonlinear response of the clamped cylindrical panel shown in 
Fig. 23 is analyzed using these models. The solutions obtained using a 4x4 grid 
of four-noded quadrilateral elements are shown in Fig. 23. The solutions obtained 
using a 2x2 grid of nine-noded quadrilateral elements are shown in Fig. 24. 

As is to be expected, the full- integration four-noded displacement model is too 
stiff. The full-integration nine-noded displacement model (with the same total 
number of degrees of freedom), though less stiff than the four-noded model, over- 
estimates the stiffness, particularly at higher loads. The mixed model with dis- 
continuous stress resultants is more accurate than the mixed models with continuous 
stress resultants developed in Refs. 17 and 18. 


ACCURACY OF FOUR - NODED QUADRILATERAL ELEMENTS 

4X4GRID 

Pq CONVERGED SOLUTION 



Figure 23 
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CLAMPED CYLINDRICAL PANEL 

ACCURACY OF NINE - NODED QUADRILATERAL ELEMENTS 

2X2GRID 

• CONVERGED SOLUTION 



Figure 24 
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USE OF REDUCTION METHODS IN CONJUNCTION 
WITH MIXED MODELS 


The use of reduction methods in conjunction with mixed models is outlined in 
Fig. 25. First, the governing finite element equations for the individual elements 
are given. Then, the vectors of fundamental unknowns (stress resultants and 
displacements) are expressed as linear combinations of a small number of vectors. 
The basis reduction and reduced system of equations are obtained in the manner 
outlined previously. It is important to note that t.hz fizdiiczd zquaid-OYU) OAZ 
quadAotic -in thz Azdaczd unknou)n& {(|)}. 


GOVERNING FINITE ELEMENT EQUATIONS FOR INDIVIDUAL ELEMENTS 



U} = VECTOR OF UNDETERMINED COEFFICIENTS 
P = LOAD PARAMETER 


REDUCED SYSTEM OF EQUATIONS 


[F] {()>} + } = p{P} 


WHERE 


[F] = E -[rjj]’’[F] [Tjj] + [Fjj] + [rjj]'^[s]'^[rjj] 

elements 


{M} 




elements 

{iji} = VECTOR OF UNDETERMINED COEFFICIENTS OR PARAMETERS 


Figure 25 
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USE OF REDUCTION METHODS IN CONJUNCTION 
WITH DISPLACEMENT MODELS 


The use of reduction methods in conjunction with displacement models, which are 
equivalent to the proposed mixed models, is outlined in Fig. 26. Note that the re- 
sulting reduced system of equations is cubic in the reduced unknowns {<|>}. The 
implication of this is that even thz. mcxzd modoJi and tkd cLL6p-tac&tmnt modet o/te 

zqul\jcitznt thoJji A.edueerf oX(l not e.qulvaZe.nt, 


GOVERNING FIN ITE ELEMENT EQUATIONS FOR INDIVIDUAL ELEMENTS 
[K] {X} + |G(X)} = p{P} 

BASIS REDUCTION 

^ ['"x] ^^^.1 • r << n 

n. r 



. a a'”* 

[•"x] = 

lx} IX} . . . . IX} 


{(Dl = VECTOR OF UNDETERMINED COEFFICIENTS 


OR PARAMETERS 
p = LOAD PARAMETER 
REDUCED SYSTEM OF EQUATIONS 
[K]|0| + {G((D)} = p IP} 

WHERE [K] = 

{G((D)} = IG (a))} 


= VECTOR OF NONLINEAR TERMS (CUBIC IN (0)) 


NOTE: EVEN IF MIXED MODEL AND DISPLACEMENT MODEL ARE EOUIVALENT, 
THEIR REDUCED SYSTEMS ARE NOT EQUIVALENT 


Figure 26 
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ACCURACY OF REDUCTION METHOD - MIXED AND DISPLACEMENT MODELS 

The nonlinear solutions obtained using the reduction method in conjunction 
with the equivalent mixed and reduced- integration displacement models are compared 
in Fig. 27 for the case of a clamped cylindrical shell subjected to uniform pressure 
loading. Seven basis vectors were generated for the unloaded shell. The variations 
of the strain energy with the loading, as predicted by the reduction method 
mixed and displacement models, are shown in Fig. 27. The high accuracy of the 
predictions of the mixed model is clearly seen in this figure. 


FULL SYSTEM (672 D.O.F.) 

▲ MIXED 
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ADVANTAGES OF MIXED MODELS OVER 
EQUIVALENT DISPLACEMENT MODELS 


The advantages of mixed models over the equivalent displacement models are 
listed in Fig. 28. These Include: 

1) S-impZ^cJjty {^ofmaZojtwn. Only quadratic and bilinear terms appear in 
the governing finite element equations. By contrast, the governing equations of 
the displacement model include cubic terms. 

2) If the geometric and material characteristics are constants within each 
element, then mo&t tk(L int(LqfiaJU> dan be. e.vaZiiated e.xadtty for the mixed elements 
(even when the element has curved faces and edges) . 

3) The mixed models are beJXeA 6Lud:ed Hie lAlith fiediLdtAjOn methodi in non- 
linear problems, in the sense that: 

a) The basis vectors are simpler to generate. 

b) The mixed models lead to higher accuracy of the solutions obtained by 

the reduced system. This is especially true for stress resultants. 


• SIMPLICITY OF FORMULATION (ONLY QUADRATIC AND BILINEAR TERMS 
APPEAR IN GOVERNING FINITE ElfMENT EQUATIONS ) 

• MOST OF THE INTEGRALS CAN BE EVALUATED EXACTLY (EVEN FOR ELEMENTS 
WITH CURVED FACES AND EDGES) 

• BETTER SUITED FOR USE WITH REDUCTION METHODS 

• BASIS VECTORS ARE SIMPLER TO GENERATE 

• BETTER APPROXIMATION PROPERTIES (HIGHER ACCURACY OF 
REDUCED SYSTEM. ESPECIALLY FOR STRESS RESULTANTS) 


Figure 28 


29 



LARGE ROTATION NONLINEAR PROBLEMS 


The third topic considered in this paper is the large rotation nonlinear 
problems. The basic features of two effective computational models are outlined 
in Fig. 29. In both models a iotdt La.QtWin.QAjan description of the deformation is 
used. Consequently, the strain-displacement relations contain tfli^onomztA^C. 
^anct^oyn o£ thz A-oteution compon&nti . 

The first computational model is a two-field mixed model with discontinuous 
stress-resultant fields at Interelement boundaries. The second model is based on 
the use of the penalty method for handling the trigonometric functions, thereby 
simplifying the analysis. 


FORMULATION 

• TOTAL LAGRANGIAN DESCRIPTION OF DEFORMATION 

• STRAIN - DISPLACEMENT RELATIONS CONTAIN TRIGONOMETRIC FUNCTIONS 
OF ROTATION COMPONENTS 

FINITE ELEMENT MODELING 

• MIXED MODELS WITH DISCONTINUOUS STRESS RESULTANTS 

• PENALTY METHOD FOR HANDLING TRIGONOMETRIC FUNCTIONS 


Figure 29 
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ELASTICA PROBLEM - FORMULATION 


As an application of the proposed computational models, consider the elastica 
problem shown in Fig. 30. In the mixed formulation, the Hellinger-Reissner two- 
field mixed variational principle is used. Transverse shear deformation, though 
small, is Included to simplify the formulation. The extenslonal strain e and the 
transverse shear strain y are trigonometric functions of the rotation (p. 

The penalty formulation, on the other hand, is based on the Euler-Bernoulll 
type beam theory with both the extenslonal and transverse shearing strains neglected. 
The axial and transverse displacements u and w are Incorporated into the functional 
through the use of constraints and penalty numbers. 


MIXED FORMULATION 


•n- = /[(Ne+Mk + Qy) 


2 2 2 
1 .n"' 4 . 4 . 1 


e, Y INCLUDE TRIGONOMETRIC FUNCTIONS OF <t> 


ic = --i 
ds 



N, M, Q ARE DISCONTINUOUS AT INTERELEMENT BOUNDARIES 


PENALTY FORMULATION 

n = / [j ^ sin(j))^ + X2 - 1 + cos iji) ds 

WHERE X2 are penalty NUMBERS 

THE RESULTING STIFFNESS MATRIX IS POSITIVE DEFINITE. 


Figure 30 


31 



ELASTICA PROBLEM - NUMERICAL RESULTS 


The displacements, rotations and total strain energy obtained by using ^ouA. 
tjjlo-nodzd oZzmdYitA and dM)o th/ie,^-node.d &tmznt6 are depicted in Figs. 31 and 32. 
Also, the deformed configurations of the beam for various values of the transverse 
load P are shown in Fig. 31. Both exact-integration displacement models (DE models) 
and mixed models with discontinuous forces (MD models) are used. As to be expected, 
the displacement models are too stiff. This is particularly true for the two-noded 
elements. By contrast, the predictions of the mixed models are highly accurate. 


LARGE-ROTATION ELASTICA PROBLEM 



— EXACT 

SOLUTION 

0 

4 DE2 ) 

DISPLACEMENT MODELS 

o 

2 DE3 ) 


+ 

4 MD2-1 ) 

MIXED MODELS 

A 

2 MD3-2 ) 




Figure 31 
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ANALYSIS AND MODELING NEEDS FOR TIRES 


The fourth part of the paper deals with future analysis and modeling needs for 
tires. The overall goal is to develop a general tire analysis capability which 
Includes (see Fig. 33): 

1) Accurate representation of the tire configuration and construction 

2) Reliable material characterization including thermo-viscoelastic response 

3) Capability for predicting the stresses and deformations due to footprint 
loading; this also includes the prediction of the contact area 

Since there is a certain degree of uncertainty in the accuracy of the various 
elements of the tire model, considerable work should be directed towards assessing 
the sensitivity of the tire response to various modeling details such as material 
characteristics, surface Inaccuracies, and variations in the tire design variables. 
The result of such sensitivity study would allow the identification of the minimum 
degree of sophistication of the model required to achieve a presclbed level of 
accuracy. 

There is also a need for identifying failure mechanisms and developing a 
verifiable failure analysis capability for tires. Use can be made of the consider- 
able experience gained in damage tolerance design concepts for fibrous composite 
structures. 

• GENERAL TIRE ANALYSIS CAPABILITY 

• ACCURATE REPRESENTATION OF TIRE CONFIGURATION AND CONSTRUCTION 

• RELIABLE MATERIAL CHARACTERIZATION INCLUDING THERMO-VISCOELASTIC RESPONSE 

• PREDICTION OF STRESSES AND DEFORMATIONS DUE TO FOOTPRINT LOADING 

(CONTACT AREA. STRESSES AND SLIP) 

• SENSITIVITY ANALYSIS - SENSITIVITY OF RESPONSE TO: 

• MATERIAL CHARACTERISTICS 

• SURFACE INACCURACIES 

• VARIATIONS IN DESIGN VARIABLES (REQUIRED FOR EVALUATION OF STRUCTURAL 

CONCEPTS AND FOR OPTIMIZATION) 

• MODELING DETAILS (IN ORDER TO DEVELOP SIMPLE TIRE MODELS) 

• FAILURE MECHANISMS AND FAILURE ANALYSIS OF TIRES 

• DAMAGE TOLERANCE 


Figure 33 
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SUMMARY 


In summary, four topics are covered In this paper; namely, recent advances in 
reduction methods for nonlinear problems, mixed models for tires, computational 
models for large-rotation nonlinear problems, and analysis and modeling needs for 
tires. (See Fig. 34.) 

Reduction methods have proven to be very effective for the nonlinear static 
analysis of structures subjected to either combined loads or displacement- 
dependent loads. However, more work is needed to realize their full potential for 
nonlinear dynamic and time-dependent problems. 

Mixed shell models with discontinuous stress resultants at element interfaces 
have high potential for nonlinear analysis of tires. These models can be easily 
Incorporated into existing general-purpose finite element programs based on the 
displacement formulation. 

Two computational models are presented for the large-rotation nonlinear 
problems. Both models use a total Lagranglan description of the deformation. 

The first model uses a mixed formulation, and the second model uses a penalty 
formulation. Both models appear to have high potential. 

As far as analysis and modeling needs are concerned, three areas have been 
identified. As is to be expected, the modeling and analysis of tires will be strongly 
Impacted by new advances in materials technology, computer hardware, software, 
integrated analysis, and CAD/CAM systems. 


• REDUCTION METHODS 

• VERIFIED FOR STATIC NONLINEAR PROBLEMS INCLUDING CASES OF COMBINED 
LOADS AND DISPLACEMENT-DEPENDENT LOADS 

• FURTHER DEVELOPMENT NEEDED FOR NONLINEAR DYNAMIC PROBLEMS 

• MIXED MODELS WITH DISCONTINUOUS STRESS RESULTANTS 

• HAVE HIGH POTENTIAL FOR ANALYZING TIRES 

• CAN BE EASILY INCORPORATED INTO EXISTING GENERAL-PURPOSE FINITE ELEMENT 
PROGRAMS 

• LARGE ROTATION NONLINEAR PROBLEMS 

• BOTH MIXED AND PENALTY FORMULATIONS PROVIDE EFFECTIVE ANALYSIS 
TECHNIQUES 

• ANALYSIS AND MODELING NEEDS 

• GENERAL ANALYSIS CAPABILITY FOR TIRES 

• SENSITIVITY ANALYSIS 

• FAILURE MECHANISMS AND FAILURE ANALYSIS 

Figure 34 
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FINITE ELEMENT METHODS FOR NONLINEAR 
ELASTOSTATIC PROBLEMS IN RUBBER ELASTICITY 


J. T. Oden, E. B. Becker, T. H. Miller, T. Endo, and E. B. Fires 
The University of Texas at Austin 


ABSTRACT 

This paper outlines a number of finite element methods for the analysis of 
nonlinear problems in rubber elasticity. Several different finite element schemes 
are discussed. These include the augmented Lagrangian method, continuation or 
incremental loading methods, and associated Rlks-type methods which have the capa- 
bility of Incorporating limit point behavior and bifurcations. Algorithms for the 
analysis of limit point behavior and bifurcations are described and the results of 
several numerical experiments are presented. In addition, a brief survey of some 
recent work on modelling contact and friction in elasticity problems is given. 
These results pertain to the use of new nonlocal and nonlinear friction laws. 


OUTLINE 

1. PRELIMINARIES AND NOTATION 

2. FINITE ELEMENT MODELS 

3. ALGORITHMS 

4. NUMERICAL EXPERIMENTS 

5. FRICTION MODELS 

6. NUMERICAL RESULTS 
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PRELIMINARIES AND NOTATIONS 


KINEMATICS 


The usual notation In finite elasticity is employed: U is the displacement 

vector, X is the position of a particle in the current configuration whose posi 
tion was X in the reference configuration, F is the deformation gradient. 
These quantities are illustrated in the figure^below. 





DISPLACEMENT VECTOR 
X - X 


DEFORMATION GRADIENT 


V(X + u) ; F„,. = 




HOMOGENEOUS, ISOTROPIC 
HYPERELASTIC MATERIAL 


We begin with a study of equilibrium problems in finite elasticity. It is 
assumed that the materials involved are hyperelastic, isotropic, and homogeneous. 
Therefore, they are characterized by a strain energy function W which is given 
as a function of invariants of the deformation tensor Q . These problems are 
complicated by the fact that nonconvex constraints must~be enforced. For compres- 
sible materials the constraint manifests itself in the condition that 0(u) = det 
F > 0 while for incompressible materials J(u) = 1. 


W = W (I-|,l2,J) = STRAIN ENERGY PER UNIT VOLUME 
I] = TRACE C 

I 2 = ^ (trace C)^ - 1 TRACE 
J = Y ^3 = det F C = f''^F 

COMPRESSIBLE MATERIALS 

J(u) = det F > 0 

INCOMPRESSIBLE MATERIALS 

J(u) = 1 
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SOME CONDITIONS ON THE FORM OF W 

Some conditions on the form of the energy W are presented below. In addi- 
tion to the fact that the energy must be form invariant under changes of the 
spatial frame of reference, other conditions must be enforced if one expects a 
well-behaved solution. Three of these are listed below; 


1. COERCIVITY; 

w(u) = W(I^(u),l2(u),0(u)) 

W(u) +°o as I |Vu| I ^ +<» 


2. SINGULAR BEHAVIOR: 


W(u) -> +<» 


1 9W I , 

IjfI 


3. QUASI CONVEXITY; 


as det F -> 0 


3 ^W(u) , , ^ n 

~ A .A .y 1 0 

9 F . a . ^ J a B 


■ta j °o 


V X , y 
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INCOMPRESSIBLE AND NEARLY INCOMPRESSIBLE MATERIALS 

« 

Two basic classes of methods are employed here to handle compressible and 
nearly incompressible materials: Lagrange-multlplier methods (mixed methods), in 
which the incompressibility constraint h(j) = 0 is accounted for using Lagrange 
multipliers and penalty methods in which the total potential energy func- 
tional n is penalized by the addition of a positive semi-definite, generally 
convex penalty functional. 


INCOMPRESSIBILITY CONSTRAINT 

h(J) = 0 

h(J) = J - 1 , - 1 , (J - 1)^, -InJ , etc. 

W = W(I^ , 12 ) - p h(J) 

p = Lagrange Mult. -Hydrostatic Pressure 


PENALTY TERMS 


n(u) 


W(u)dX - f(u,v) 


a 


+ E 




g(J(u))dX 


g(J) > 0 g(J) = 0 J = 1 
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EQUILIBRIUM CONDITIONS 


The finite elements employed are based on various alternative statements of 
equilibrium conditions for elastic bodies: 


1) ENERGY FORMULATION 


n(u> = TOTAL POTENTIAL ENERGY 


= W(u)dX + f(u,v) 

Jfl - - ~ 

f(u»v) = f f(u)*vdX + f t(u)*vds 


n(u) < n(v) for all V in 
K = {v : I f W(v)dX|<* ; 

V = 0 on ; det 7v > 0 or = 1} 


2) AUGUMENTED LAGRANGE 


L(u,F,X) = n(u) + I [ |Vu - Fl^dX 
- = = - = 

+ f X ; (ffu - F)dX 

>a- ' ' 


L : V X {G: det G = 1 ) x { ■ ■ )* 


The augmented Lagrange method combines features of both penalty and 
Lagrange-multiplier schemes. The Incompressibility constraint can be satisfied 
a priori in a straightforward manner for Mooney-Rlvlin materials, with the 
result that the method is extremely fast and efficient when used in conjunction 
with, e.g. , Uzawa's method (ref 1). (See equation.) 


3. VIRTUAL WORK I 


! aw(u) r 3h 30 

- J„“3J' ■ ''i.ci'* 

= f(uiV) for all V in V 
|jjqh(J(u))dX . 0 for all q in Q 


VIRTUAL WORK II 


f V. dx f o-' f 

'a 3u. >a 

■Lfd 


dg dO 

7J * 3u. 

-c,a 


f(u ,X) for all V in V 


-1 3g(J(u^)) 
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FINITE ELEMENT METHODS 


SOME ELEMENT FAMILIES 

The figure below Illustrates some of the standard finite element methods 
employed. Unfortunately, not all of these methods are numerically stable since some 
may not satisfy the generalized LBB condition of Oden and LeTallec (refs. 2 and 3). 
This condition is given in the equation on the following page, where Q is the space 
of Lagrange multipliers, Vj^ is the space of finite element approximations of dis- 
placement, p is the hydrostatic pressure and | IVv^^l |o,p is an appropriate energy 
norm on y- . The existence of a 3h ^ 0 lndependent~of mesh size h is necessary for 
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THE GENERALIZED LBB CONDITION 
(refs. 2 and 3 ) 



The numerical stability of mixed- and penalty finite element 
methods depends upon this condition, and particularly on the behavior 
of the parameter vith the mesh size h. 
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STABILITY RESULTS 


Oden and Jacquotte (ref, 4) have recently completed a mathematical analysis of 
the LBB condition for incompressible viscous flows. Some of these results appear to 
be directly applicable to the finite elasticity problem. The following table summer 
Izes the behavior of certain finite element methods for these constrained problems. 
The results generally fall into three categories: the stability parameter 3h is 

independent of h and the method is stable, Bh is dependent on h and the method 
is unstable having spurious pressure modes, or the element is "locked", meaning that 
the penalty parameter £ depends on the mesh size and that the displacements 
approach 0 as e tends to zero for a fixed mesh size h . 


VELOCITY APPROX. 


QUADRATURE RULE 
(PRESSURE APPROX. Q ) 


CONVERGENCE 

RATE 







































ALGORITHMS 


Several different algorithms are employed for the analysis of finite elasti- 
city problems. These include the augumented Lagrange/Uzawa methods, continuation 
(incremental loading) methods, and homotopy methods. In the present paper, augmented 
Lagrange methods are discussed briefly, but the focus is on the continuation- type 
methods of, for example, Rlks (ref. 5), Keller (ref. 6), Crisfleld (ref. 7), and 
Padovan (ref. 8). 

1. AUGUP1ENTED LAGRANGE/UZAWA 

1. EXTREMELY SIMPLE & FAST 

2. FOLLOWS STABLE BRANCHES 

3. NOT EASILY ADAPTED TO GENERAL MAT'LS 

2. CONTINUATION (INCREMENTAL LD'G) 

1. Kiks (ref. 5)j Wempner (ref. 9) 

2. Keller (ref. 6).. Rheinbolt (ref. 10) 

3. Crisfield (ref. 7)j Padovan (ref. 8) 

3. HOMOTOPY METHODS 

f(x,p) = 0 

9f (x(s ) ,p(s ) ) . af (x(s ) ,p(s ) ) . 

2 X + — — p = 0 

ax ~ ap 

N(x(s),p(s)) = 0 

+ ODE SOLVER 
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AUGUMENTED LAGRANGE METHOD 
The augumented Lagrange is a super-fast method. 


AUGUMENTED LAGRANGE BLOCK RELAXATION 

1 . - P (Vu" - f") 


2 . 


L(u" , . X") < L(u[J . G . X") 

ISOLVED EXPLICITLY, N + 2,3| 


^v^^~k+1 ’ ~k+l ’ ~ ^ ^ 


CONTINUATION METHODS 


f(x(s) , p(s)) = 0 

k{s) . x(s) + p^(s) = 1^ N(x,p) = 0 
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NUMERICAL EXAMPLES 


Consider the inflation of a thick rubber spherical shell subjected to external 
loading. The shell is pressurized until it snaps through. Pressure continues to 
increase until the shell inflates. The material is assumed to be an incompressible 
Mooney-Rivlin material. The numerically stable Q2/P1 element is employed. Good 
results are obtained using the Rlks-Crlsfleld method (refs. 5 and 7) with Newton- 
Raphson correction. Geometry of the shell is shown below and this is followed by 
several figures which illustrate numerical results. It is noted that the stiffness 
of the shell is strongly dependent on the material properties. In particular, for 
a fixed value of the Mooney-Rivlin constant Cj, if C 2 is chosen sufficiently small, 
a limit point type behavior occurs in the inflated shell. This represents the phe- 
nomenon of a large decrease in pressure in an inflated balloon with a large increase 
in strain at the crown. This limit point behavior disappears for larger values of 
C2- (See ref. 11.) 
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CALCULATION OF BIFURCATION PATHS 


The Riks method (ref. 5) makes it possible to calculate bifurcation paths as 
well. A second example Involves the buckling of a thick rubber plate under end thrust. 
This problem was analyzed by Sawyers and Rivlin (ref. 12), and provides a good example 
to test the bifurcation capabilities of the code. For low initial aspect ratios, a 
barreling mode of bifurcation is obtained, whereas for thinner slabs a flexural mode 

is obtained. The following figures show computed numerical results for these ele- 
ments. Bifurcation branches were calculated for this problem as well. The fol- 
lowing figures show the results of these calculations. After checking the energy 
of the system on each branch, the branch of lower energy is then calculated . In 
the first figure, branch one has lowest energy and branch two indicates a second- 
ary bifurcation. The lower branches in each of the subsequent figures also rep- 
resent lower energy equilibrium states. (See ref. 13.) 


PLATE UNDER END THRUST 
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first three buckling modes for a 


PLATE UNDER END THRUST 




A = 0.892 

X 
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BUCKLED CONFIGURATIONS OF A PLATE FOR FIRST MODE BUCKLING 


ORIQ9ML PMGE IS 
$fi mR QUAUTY 



AL/L =0.30 
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COMPUTATIONAL PROBLEMS DURING HIGH COMPRESSION 


The augumented Lagrange method failed for problems with high compression as 
indicated in the figure below. An oscillatory mode was obtained at compressions of 
30%. Here the Ql/PO element was employed. In view of this, the continuation tech- 
niques with Q2/P1 elements were attempted in the following examples. 


30% COMPRESSION 
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FIXED END CYLINDER 


Area shown in 
defornatlon plots 




- N 











COMPRESSION OF A RUBBER BLOCK 

Here a rubber block of revolution is one-quarter of the cross section of a 
block of revolution is shown. The block is subject to compression by uniformed dis- 
placement of the edges. A strong singularity is developed in the free edge leading 
to a cusp as indicated in the calculated deformed configurations. Ultimately, contact 
of these surrounding regions occurs. This is a very difficult class of elastostatlc 
problems and contact conditions must be incorporated in the analysis procedure to 
handle these problems. 


PLANE STRAIN BLOCK 
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enlarged view of deformed configuration 



At limit point 



After limit point 



Element boundary overlap at the cusp 
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NEW FRICTION MODELS 


It is now the generally held view that the classical Mohr-Coloumb law of friction 
is inadequate on both physical and mathematical grounds when modelling friction 
effects in real materials. For this reason, a preliminary study of several new non- 
classical friction laws have been undertaken. These friction laws Include nonlocal 
effects to approximate deformed asperities, adhesion effects, and elastic junction 
effects, and ultimately will Include fracture and damage models of the junctions on 
the contact surface. The general form of these contact laws is given below, where 
aj is the tangential frictional stress, V is the coefficient of friction, which 
may depend upon the number of loading cycles and the stress state on the contact sur- 
face, gradients of deformation, etc., Sp is a smoothing operator with p represent- 
ing a characteristic dimension of deformed asperities on the contact surface, 
is the normal contact pressure, and 4>g. represents a compliance function which models 
the elasticity and elastoplasticity of interface junctions. This compliance reaction 
is given as an anti-symmetric function of the value of the tangential velocity vec- 
tor Uj on the contact surface. An algorithm has been developed for implementing this 
new friction law. Some preliminary results are indicated in the final two figures 
of this paper. (See also refs. 14 through 16.) i 



a.|.(u) = v(Vu,t) S (aj^(u)) ( I I Uj | 1 ) 


! I 


CYCLIC/DAMAGE NONLOCAL JUNCTION ELASTICITY & PLASTICITY 
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NEW FRICTION MODEL APPLICATIONS 

^ example of the indentation of a rigid cylindrical stamp into an elastic 
slab is^considered. A non-classical friction law is used, with p = -0.1, 

^ ~ 1000, y = 0.6. Nine-node bi-quadratic elements were used. Deformed 
shape and stresses on the contact surface are shown. 
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NEW FRICTION MODEL APPLICATIONS (concluded) 


R 


1000 J 


750 J 


tl 


Normal Contact Pressure 
“A- Tangential Stress 


500 






250 


o- — o-*‘ 


11 


A ' 
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T 
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Vi 
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APPLICATION OF TIRE DYNAMICS TO 
AIRCRAFT LANDING GEAR DESIGN ANALYSIS 


Raymond J. Black 

Bendix Aircraft Brake and Strut Division 


ABSTRACT 


The tire plays a key part in many analyses used for design of aircraft landing 
gear. Examples include structural design of wheels, landing gear shimmy, brake 
whirl, chatter and squeal, complex combination of chatter and shimmy on main 
landing gear (MLG) systems, anti-skid performance, gear walk, and rough terrain 
loads and performance. This paper discusses which tire parameters are needed 
in the various analyses. 

Two tire models are discussed for shimmy analysis, the modified Moreland 
approach (ref. 1) and the Von Schl ippe-Dietrich approach (ref. 2). It is shown 
that the Moreland model can be derived from the Von Schlippe-Dietrich model by 
certain approximations. 

The remaining analysis areas are discussed in general terms and the tire 
parameters needed for each are identified. The conclusion of the paper is that 
accurate tire data allows more accurate design analysis and the correct pre- 
diction of dynamic performance of aircraft landing gear. 
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INTRODUCTION 


A number of important design analysis areas of aircraft landing gear require 
accurate data on tire dynamic parameters. One of the more obvious examples is wheel 
design. The primary loads acting on the wheel act through the tire onto the wheel 
bead seat area. The wheel cannot be efficiently designed without accurate data on 
tire to wheel contact loads for the entire spectrum of operating conditions to which 
the landing gear system is subjected. This includes landing impact, braked roll, 
taxiing, turning, cross wind landing, and in some cases obstacle impact. The latter 
is particularly true for carrier based airplanes. The wheels of such aircraft must 
accommodate running over arresting cables when the tire is in the near bottomed con- 
dition. Modern aircraft wheel design makes use of experimentally and analytically 
determined bead contact forces to arrive at optimized designs of the wheel cross 
sections. In spite of increased demands for more durable wheels having longer 
fatigue life and increased capabilities, wheel weight per unit load has decreased in 
recent years. 

Other design analysis areas involving the entire landing gear system include: 

• Landing gear shimmy (both nose and main landing gear) 

• Brake-landing gear whirl stability 

• Brake chatter and squeal 

• Complex combinations of chatter and shimmy 

• Anti -skid performance 

• Gear v/al k 

• Rough terrain operations 

This paper will concentrate primarily on the tire parameters needed in the first four 
of the above problem areas. 


72 



LANDING GEAR SHIMMY 


There are two tire models commonly used for aircraft landing gear shimmy 
analysis: the modified Moreland model (refs. 1, 3, and 4) and the Von Schlippe- 
Dietrich model (refs. 2, 5, and 6). The equations for the forces and moments at 
the tire footprint are similar for the currently used versions of the two models and 
as a result there are a large number of parameters which are common to both models. 
Both models require that the vertical load deflection characteristics of the tire be 
known along with the general load, inflation pressure and speed envelope of opera- 
tion of the aircraft. In addition the following parameters must be known for the 
operational envelope: 

m^ = mass of tire (kg) 

2 

Ipt ~ polar moment of inertia of inflated tire (kg*m ) 

2 

I^t “ diametral moment of inertia of inflated tire (kg-m ) 

= lateral spring rate of tire (N/m) 

C|^ = lateral equivalent viscous damping coefficient (N-sec/m) 

l-L^ = rolling tire torsional spring rate (N-m/m) 

I_Iq = rolling tire torsional damping coefficient (N-m/sec) 

Rg = loaded radius of the tire (m) 

= slope of drag force versus slip ratio for small slip ratios (N) 

R^ = tire rolling radius (m) 

For the modified Moreland model, two additional parameters are required: 

C = coefficient of yaw (rad/N) 

= yaw time constant (sec) 

These parameters are used in Moreland's cornering relationship: 


CF 


nj 


= ^ 


tj 






( 1 ) 
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where 


j = index indicating the jth wheel of a multi-wheel landing gear 
F . = lateral force acting on tire footprint (N) 

1 1 J 

= yaw angle or sideslip angle of tire footprint with respect to 
^ the tire (rad) 

Dots over a variable are used to indicate differentiation with respect to time. 

In the Von Schl ippe-Dietrich model it is shown that the cornering characteristics 
come about due to continuity between the tire distortion in the footprint's 
most forward position and the lateral distortion of the tire as a function of 
peripheral distance, combined with adhesion of the tire footprint centerline to the 
ground surface. Denoting the lateral displacement of the tire outside the foot- 
print as X and the peripheral distance around the tire as s (s = 0 corresponds to 
the forward tip of the footprint and s > 0 corresponds to moving upward around the 
tire center) it is found from experiments that when a pure lateral load is applied 
to the tire X is given by 


X = 



( 2 ) 


where Xp = lateral deflection of footprint (constant within footprint) (m) 
s = peripheral distance (m) 

a^ = experimental constant between zero and 1 (dimensionless) 

L = experimentally determined "relaxation" length (m) 

For type VII tires, aj is very close to 1. However, for many new design tires 
with a relatively small aspect ratio and a relatively small width between wheel 
flanges, a^ is significantly less than unity. For example, on a H45 x 17-20, PR22 
tire, typical values of a^ range from 0.80 to 0.92, depending on the vertical load 
and the tire pressure. 

As a consequence of Equation (2), it can be shown that when the wheel moves forward 
an amount ds the most forward position of the tire footprint will move an amount 


dX 


1 




(3) 
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relative to its position prior to rolling forward. The quantity a^^/L is denoted by 

c = 3,-^/Jj (4) 

It can be shown that as a consequence of Equation (2) the steady state "cornering 
power" or cornering coefficient is given by 

N = % ^ (l + j (5) 

where 

N = cornering coefficient (M/rad) 

= tire coefficient defined by Equation (4), (1/m) 

h = half length of the tire footprint (m) 

For transient conditions the general kinematic situation is illustrated in Figure 1. 
Due to strut deflection, the wheel is tilted at an angle and twisted in the yavif 
direction at an angle xj/'s , The coordinate of the wheel axle center is ya and the 
coordinate of the most forward position of the tire footprint is y^p* The wheel is 
rolled forward an amount ds during which time the position variables of the tire 
change by amounts d(bL, d i// ' dy,, dX , and C^Xds as defined by Equations (3) and 

/ A \ D b ^ 



Figure 1.- Tire distortion nomenclature. 
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With this in mind we can, from Figure 1, write the quantity + dy^„ in two 
ways. We first consider the path a to b to c to d to e which gives^ ^ 


7fp + dy^p = Ya + Rg<^B + dy^ + Rgd0' - h(;^; + - (X + dX) (6) 


Next, we consider the path a to b to f to g to i to j which gives 


+ Rg0B - h 03 ' - X - + C^Xds (7) 


Equating (6) to (7) gives 

dy^ + Rgd0g - hdl//g'- dx = - %'ds + C^Xds 


( 8 ) 


which can be written as 

dX 


^ + C.X = T-^ + Rg — ^ + 0s h ^ 

ds X ds S ds ® 


(9) 


We 
from 


can eliminate X from Equation (9) by defining X in terms of y^ and y^ . Again 
m Figure 1, 

X = y^ - yfp + Rg0' - h0/ 


Substituting (10) into (9) gives 


dyfp 

~dT +^xyfp 


Cx(ya + Rg0s) - <1 + Cxh)0s' 


Or, letting -~ = D, this substitution will yield 
ds 


fp D +c 




- ( 


1 + C^h 


)^s'l 


( 11 ) 


( 12 ) 


As a result of adhesion of the tire footprint to the ground, we can write 



Taking the series expansion for the left hand side of Equation (13) gives 


yfp(s = -h) = yfp - h 


ds 


,2 


l£p 

2 1 ds^ 


(14) 


1 

Applying (14) to (12) and making the substitution D = — where = d/dt gives 


/ n »» 0 


3!v 


3-0^3 + . . .)(cx(y, + Vb> - + “ 

+'=x)k + ><»<#> B -^) 


(15) 


Equation (15) is the basic cornering relationship used in the current version of the 
Von-Schlippe-Dietrich tire model. 


Although Equation (15) does not appear to resemble the Moreland cornering relation- 
ship, it can be shown that they are similar. If the infinite series representing 
the distance delay is restricted to the first two terms, 1 - (h/v)D , and the fol- 
lowing kinematic relationship is introduced: ^ 


= -®t(ya + Rg0B - 


(16) 


where V is the forward velocity of the airplane, then Equation (15) can be 
rewritten as: 


/I + hc>,\ u / . (i + hc^\ . \ 


(17) 


Since the side force is given by 


+ c^A 


(18) 


Equation (17) can be written 


K^(l + hC;^ ) 


Ct • C,. h(Ki + CTDf) 

F„ - +-L — I — lA- 

•'i 


VKj^(l + hC^ 


t) / • 


(1 + Cy^h) 


(19) 


For large values of V the far right hand term is small. Thus Equation (19) is 
similar to Equation (1) with 


( 1 + hC^) 


~ C or — from Equation (5) 
N 


(20) 


and 


( 21 ) 
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From the latter it is apparent that Moreland's cornering relationship is equiv- 
alent to an approximation of the Von Schlippe-Dietrich model. This explains why the 
two models give similar results when used in shimmy analysis. It has been found 
that in many cases the Moreland model is slightly more unstable than the Von 
Schlippe-Dietrich model. It is often more conservative to use the Moreland model 
for landing gear design. Moreland's model is also more simple from a computational 
standpoint since one does not have to deal with the infinite series or, equivalently, 
the time delay function represented by Equations (14) and (13) respectively. 

The forces acting at the tire footprint are only part of the shimmy problem 
formulation. Figure 2 illustrates the factors considered in the total problem, 
namely the airplane, the local structure between the airplane and the landing gear, 
the landing gear structure, and finally the footprint of the tire. 

The fuselage is represented by the three dimensional motion of a fuselage 
reference point. The three motions are a lateral motion yf, a roll motion 4>f, and a 
yaw motion;//^. Three input forces corresponding to these three directions act from 
the landing gear. The motion of the fuselage is given by solution to model equa- 
tions covering the first several normal modes of the fuselage. The equation of 
motion of the fuselage is: 


[m] j*qj + (1 + i 77 f ) [k] jqj = [ 5 ] 


T \]ya 
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where 

[M] = fuselage modal mass matrix (kg) 

[Kl = fuselage stiffness matrix (N/m) 

= loss factor for fuselage structure 

[6] = modal deflection matrix 

forces and moments in the y^xf/and 0 directions 

q = vector of n normal modes 



To obtain the fuselage deflection at the reference point the q vector is multiplied 
by the transpose of the modal deflection matrix. 



The relative deflections due to the structure from the fuselage reference point to 
the landing gear attach point are given by: 
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where 


^t-^f 




r= relative deflections in the y, i// and 0 directions 


[f] = flexibility matrix 


y 

V 

T0 


forces and moments in the y, \}/ , and 0 directions 


The mathematical model of the landing gear is fairly straightforward. Care must 
be taken in properly accounting for cross coupling of flexibilities and inertial 
cross coupling brought on by off-center masses of the landing gear system such as 
torque arms, steering cylinders, lights and attachment fittings. Since each gear 
differs in details the equations of motion will not be given here. 

The remaining equations are those of the tire forces and kinematic restraints. 
For both the Moreland and Von Schlippe-Dietrich models, these equations are: 


F^. = K A . + n A . (25) 

nj 1 j L j 

for the side force, 

^ (26) 

for the moment, and 

V( 0' - = -y^ - Rg0g + Aj (27) 

for the kinematic constraint. 
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For the Moreland tire model the tire equations are completed using Equation (1) as 
the final relationship while for the Von Schlippe-Dietrich model Equation (15) is 
used. 


For dual wheel landing gear systems there is an additional moment acting from 
the ground about an axis normal to the ground plane. This moment, M^^, is governed 
by the following differential equation: 


M + (c R R„/VI ) M + 
w w r g p w 


2K /I /m ( 

w P / w 


dt 


(C„l//2V) - K„C„l//VI (28) 


where in addition to previously defined parameters 

2 

Ip = polar moment of inertia of one wheel and tire (kg-m ) 

Ig = distance between v/heel centerlines (m) 

= torsional spring rate of corotational axle (N-m) 

If = 0 then Equation (28) still may be used by simply setting = 0 and retain- 
ing the remaining terms. It can be noted that an additional tire parameter, is 
needed to deal with dual wheel landing gear systems. This parameter defines the tire 
drag forces in terms of the percent tire slip for small slip ratios. 

Figure 3 shows the Space Shuttle nose landing gear (NLG), which was analyzed 
extensively for shimmy stability using computerized solutions of the equations of 
motion previously described, dynamometer testing, and "runway" tests on the Langley 
Landing Loads Track. This work was extremely important since there were no high- 
speed taxi tests in the Space Shuttle development program. The first "test" of 
the landing gear was an actual landing at over 200 mph. 

The computer study was used to select the optimum damping for the steer damp 
unit, to determine sensitivity to wear, friction, and tire parameters, and generally 
to establish what the margin of stability was for the landing gear system. This was 
then verified and improved upon by the dynamometer test. Application to a flat 
runway was verified by comparison of dynamometer and Langley Landing Loads Track 
test results. The Langley Landing Loads Track could not be used for the full range 
of speeds and vertical loads because of limitations on track capacity. The dyna- 
mometer tests covered the full range of speeds and loads. 

Figure 4 demonstrates, for a typical case, the correlation between the computer 
analysis and the dynamometer (120" flywheel) results. Figure 5 shows the correla- 
tion between the dynamometer tests and the NASA dry concrete runway data from the 
Langley Landing Loads Track. Both sets of results correlate closely, showing that 
the analytical model can be used with confidence for prediction of the Space Shuttle 
NLG's stability over a wide range of landing and roll-out conditions. 

The process of test and analysis for shimmy stability is described in Reference 

7. 
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MAIN LANDING GEAR SHIMMY 


In addition to NLG's, shimmy can be encountered on main landing gear systems 
(MLG's). Three very popular commercial jet aircraft having dual wheel main landing 
gears experienced similar MLG shimmy problems when they were first introduced. The 
problem was traced to a side brace configuration which caused twist in the landing 
gear due to application of a side load. The side braces were all connected to the 
forward part of the shock strut for the three different aircraft. Modifications 
were made to the MLG's to stabilize the systems. In two cases a torque arm damper 
was introduced while in the third case a negative (forward) mechanical trail was 
used to stabilize the shimmy. 

In general MLG's are more complex than NLG's; however, the same tire para- 
meters are needed for MLG stability studies as are used for NLG's. The complexity 
is associated with the landing gear structure itself. This is illustrated in 
Figure 6 which is a drawing of the six wheeled C5A MLG. The analysis of this system 
used a 48 degree of freedom analysis reduced to a twelve degree of freedom system 
for non-linear simulation studies. The 48 degrees of freedom were necessary to 
develop an accurate mass and stiffness matrix for the reduced system. 

Sample stability results are given in Figure 7 and illustrate the importance of 
accurate tire parameters. Here the only difference between the two runs is the yaw 
time constant. The solid curve corresponds to a yaw time constant (C^^) of 0.06 
seconds while the dotted curve applies to a Cj value of 0.03 seconds. The actual 
yaw time constant for the 49 x 17 tires is considerably lower than both of these 
values but the comparison served to show the sensitivity of shimmy stability to tire 
parameters. 



Figure 6.- C5A main landing gear. 
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The net result of the shimmy stability study showed the C5A MLG to be stable, a 
result which was later verified by in-service experience. The analysis was impor- 
tant since this was the first time that this six wheeled landing gear configuration 
was used on an aircraft. 

The diversity of MLG designs is illustrated by the complex single wheel MLG 
shown in Figure 8. This was a preliminary design for the gear which was complicated 
by the fact that it had to retract around a wing mounted rocket. Needless to say no 
single set of equations of motion suffice to cover all landing gear systems. 
Fortunately, however, the tire parameters previously given apply to all landing 
gears. 



'□on .K .10 .IS .» .JO .K .no .15 .iO 

Time, Second 


Figure 7.- Comparison of transient motion for two values 
of the yaw time constant (C5A MLG). 



Figure 8.- Single wheel MLG design for a fighter aircraft. 
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BRAKE WHIRL 


Brake whirl is a highly destructive self excited vibration of the wheel, tire, 
brake and landing gear structure. It is caused by a localized high pressure region 
in the brake stack that brings about out of plane bending of the axle. The result- 
ing motion is characterized by a sinusoidal bending of the axle which takes place in 
two planes, the z-y plane and the x-y plane, as shown in Figure 9. The two motions 
are 90° out of phase which makes it appear that the bending motion is whirling about 
the axle centerline at a speed corresponding to the radial natural frequency. 

The mathematical analysis of the brake whirl problem requires that a number of 
tire properties be known. Since the frequency of the whirl motion is generally in 
the range of 200 to 400 Hz for most landing gears the tire does not respond as it 
does for low frequency motions. The tire tread remains fixed in the x-z plane except 
for its steady forward motion and rotation. The heavy bead area of the tire moves 
with the wheel adding to its effective mass and its moments of inertia. The region 
between the tread and the wheel is treated as a spring with radial and angular 
spring constants and associated loss factors. 


y 


Figure 9.- Brake whirl. 
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The specific tire parameters are: 


K = spring rate of tire associated with radial displacement of the wheel 
with the tread held fixed (N/m) 

7]^ =. loss factor associated with 

= diametral angular spring rate of the tire associated with twist of the 
wheel with the tread held fixed (N*m/rad) 

\r ~ factor associated with 

Mtw = tnass of the tire moving with the wheel (the bead area of the tire) (kg) 
Ip = polar moment of inertia of m^^ (kg-m ) 

'2 

Idt = diametral moment of inertia of (kg*m ) 

= tire rolling radius (m) 


Brake whirl stability is a somewhat more serious problem on carbon brake systems 
because of their mass distribution properties. A careful analysis as well as design 
changes to stabilize the motion are important in the initial stages of the brake 
development. 

Measurement of the tire stiffness parameters is difficult. As a first approx- 
imation is taken to be the average of the fore and aft tire stiffness and the 
vertical stiffness of the tire. The value of is calculated based on lateral de- 
formation of the tire at both the footprint and 180° from the footprint. The loss 
factor is estimated to be approximately 0.15 for typical aircraft tires. The mass 
Mtw ''s roughly 30% of the tire for typical aircraft tires. Calculation of Ip and 
I(jt can be carried out knowing M-^w and the location of the bead bundle. The^rolling 
radius can be determined by conventional methods and is equal to the undeflected 
radius of the tire minus one-third of the tire deflection as a very close approxi- 
mation. 
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BRAKE CHATTER AND SQUEAL 


Brake chatter and squeal are self excited motions of the brake and landing gear 
system brought on by so called negative damping at the rotor-stator interface. The 
negative damping results from an increase in brake torque for a decrease in the in- 
stantaneous brake slip velocity. This characteristic yields a term in the equations 
of motion which looks like a damping term in that it is a force or torque propor- 
tional to velocity, but has a sign opposite to that of the dissipative damping of the 
motion. Thus, the lining characteristic has been termed negative damping. 

The chatter motion is a low frequency motion (5 Hz to 15 Hz) and consists of 
fore and aft bending motion of the shock strut. Squeal consists of wind-up of the 
stationary parts against the torque take-out system and has a natural frequency of 
150 to 500 Hz. 


A relatively simple model of brake chatter and squeal uses three degrees of 
freedom. These are the fore and aft bending motion of the landing gear shock strut, 
the torsional wind-up of the nonrotating parts of the brake against the torque take- 
out system, and the forward velocity of the airplane. The nonrotating parts of the 
brake consist of the brake piston housing, the pressure plate, the backing plate, 
the stators of the brake, and the torque tube of the brake which carries the 
stators. The torque take-out system consist of an axle flange on most single and 
dual wheel landing gears and a brake equalizer rod on truck type landing gears com- 
monly used on larger aircraft. 


The tire parameters used in the simplified chatter and squeal analysis consist 
of the following: 

R^ = tire rolling radius (m) 

Rg = loaded radius of tire (m) 

= slope of drag force versus slip ratio for small slip ratios (N) 

= mass of tire (kg) 

2 

Ip = polar moment of inertia of tire (kg-m ) 

C = ratio between the shift in the center of pressure and the aft deflection 
^ of the footprint 


This three degree of freedom analysis is quite simplified but is useful for 
preliminary studies of stability of chatter and squeal, particularly for the 
selection of lining materials that will not result in self-excitation of the vi- 
bratory motions. At a later stage in the landing gear development, a more detailed 
analysis of brake chatter and squeal such as the one that is described in the next 
section can be carried out. 
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COMPLEX COMBINATIONS OF CHATTER AND SHIMMY 


For most modern main landing gear systems the chatter motion is not simply a 
fore and aft bending action of the shock strut but rather consists of a complex 
motion of the gear involving fore and aft, lateral, roll, and yaw motions of the gear. 
Thus many of the same tire parameters that are used in shimmy analysis are also 
needed for an analysis of the brake chatter and squeal. 

Figure 10 shows a typical truck type main landing gear system and some of the 
nomenclature and motions considered. There are ten degrees of freedom in the 
model : 

Xg = aft deflection at the truck pivot point (m) 

Vg = lateral deflection at the truck pivot point (m) 

0g = roll angle of truck (rad) 

0 g = end bending angle of shock strut (rad) 

0g = yaw angle of truck (rad) 

= angular velocity of one of the forward wheels (rad/sec) 
cOg = angular velocity of one of the aft wheels (rad/sec) 

^sl ~ angular rotation of the forward brake stationary parts (rad) 

©ss = angular rotation of the aft brake stationary parts (rad) 

Xp = forward translational motion of airplane (m) 

A fairly extensive set of tire parameters is needed for the model. The shimmy 
action uses the simple cornering relationship equating the tire lateral force to a 
sideslip angle divided by the coefficient of yaw. This is equivalent to the modified 
Moreland model with C^ = 0. 
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The tire parameters used in the model are: 

= fore and aft tire spring rate (N/m) 

= vertical tire spring rate CN/m) 

C = coefficient of yaw (rad/N) 

lifP = lateral coefficient of friction between tire and runway, a function 
of slip ratio 

^\f/m ~ niaximum value of 

11 ^ = ratio between self-aligning torque and sideslip angle (N-m/rad) 

= maximum fore and aft coefficient of friction between tire and runway 

C = ratio between the drag force and the slip ratio for small slip 
ratios (N) 

= ratio between the shift in the center of pressure and the aft 
deflection of the tire 

= mass of tire (kg) 

2 

= polar moment of inertia of inflated tire (kg>m ) 

2 

Idt = diametral moment of inertia (kg-m ) 

Rg = loaded radius of tire (m) 

= rolling radius of tire (m) 


A sample result of the nonlinear analysis is shown in Figure 11. In this run 
the negative damping from the brake lining was set relatively high and the system 
damping, made up of several nonlinear terms, was set at about one half the expected 
value. A pulse was added to the torque to cause excitation at the midpoint of the 
10 mph stop. It can be seen in Figure 11 that the initial motion is stable but 
that the pulse causes enough disturbance that the resulting motion is near neutral 
stability. The mode shape of the motion is far from a simple fore and aft truck 
motion. Note that the lateral motion of the truck is the same order of magnitude as 
the aft motion of the truck. Also, there is some yaw motion and considerable roll 
motion of the truck. 

Studies of this type help to establish the amount of negative damping that can 
be tolerated by both the chatter motion and the high frequency squeal motion of the 
landing gear. Also mode shapes and frequencies of the several low frequency chatter 
motions can be established together with the margin of stability of each mode. 
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Figure 10.- Nomenclature and motions for complex chatter 
and shimmy analytical model. 
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Figure 11.- Sample computer run for MLG chatter and shimmy analytical model 



OTHER ANALYSIS AREAS 


The other analysis areas mentioned at the beginning of this paper also require 
detailed information on tire characteristics. Anti-skid performance requires a 
model of the brake pressure to brake torque transfer function along with the 
dynamic characteristics of the anti -skid system. The remainder of the model is 
similar to that used for chatter and squeal. Gear walk which is a stick-slip 
motion similar to chatter can come about due to stick-slip action at either the 
ground surface or the brake interfaces. Rough terrain operations require a tire 
model capable of predicting vertical and drag loads due to obstacle or hole impact 
at obstacle lengths which range from a small fraction of the footprint length up to 
several times the footprint length. The short lengths present more of a problem but 
methods are available for predictions of tire loads for all lengths. Rough terrain 
operations analysis also requires a detailed model of the shock strut stroking 
dynamics in addition to the tire model. 


CONCLUSIONS 


Tire dynamic parameters play a key part in many design analysis areas for air- 
craft landing gear systems. Advanced analytical methods to predict tire parameters, 
such as finite element methods, are needed together with experimental methods for 
verification of predictions for selected cases and determination of parameters that 
cannot be determined analytically. Accurate tire data allows more accurate design 
analyses and correct prediction of dynamic performance of aircraft landing gear. 

The net result is a more reliable and more efficient system. 
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A TIRE CONTACT SOLUTION TECHNIQUE 


John T. Tielking 
Civil Engineering Department 
Texas A&M University 


EXPANDED ABSTRACT 

An efficient method for calculating the contact boundary and interfacial pres- 
sure distribution has been developed. This solution technique utilizes the discrete 
Fourier transform to establish an influence coefficient matrix for the portion of 
the pressurized tire surface that may be in the contact region. This matrix is 
used in a linear algebra algorithm to determine the contact boundary and the array 
of forces within the boundary that are necessary to hold the tire in equilibrium 
against a specified contact surface. The algorithm also determines the normal and 
tangential displacements of those points on the tire surface that are included in 
the Influence coefficient matrix. Displacements within and outside the contact 
region are calculated. 

The solution technique is implemented here with a finite-element tire model 
that is based on orthotropic, nonlinear shell of revolution elements which can 
respond to nonaxls 5 nnmetric loads (refs. 1, 2). The basic characteristics of this 
relatively comprehensive tire model are described in reference 3. This presentation 
will focus on the contact solution technique published in reference 4. A sample 
contact solution is presented for the 32 x 8.8 Type VII aircraft tire that was 
studied in reference 5, 
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FINITE-ELEMENT TIRE MODEL 


The tire is modeled by an assemblage of axisymmetric curved shell elements. 

The elements are connected to form a meridian of arbitrary curvature and are located 
at the carcass midsurface. Figure 1 shows the assembly of 21 elements along the 
midsurface of a G78-14 tire, for which calculated results are shown in this paper. 

A cylindrical coordinate system is used, with r, 0, and z indicating the radial, 
circumferential, and axial directions, respectively. Each element forms a complete 
ring which is initially axisymmetric with respect to z. The elements are connected 
at nodal circles, hereafter referred to as nodes. 

The finite elements are homogeneous ortho tropic with a set of moduli specified 
for each individual element. The orthotropic moduli for each element are determined 
by the ply structure surrounding the element. Each ply (on each element) is speci- 
fied separately, thereby allowing the model to Include carcass details such as an 
overhead belt, sidewall reinforcement, and turnups. A turnup is Included in the 
G78-14 tire model. It was found necessary to Include the turnup in the model to 
obtain the correct inflated shape. 
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SINGLE HARMONIC RING LOADS 

The finite element tire model will respond to single harmonic ring loads on 
the nodal circles. An approximately linear load-deflection response is obtained 
when an Individual ring load is applied to any node of the pressurized tire model. 
An example ring load-deflection calculation for the G78-14 tire model is shown 
in figure 2. A harmonic sequence of stiffness matrices is obtained by applying a 
sequence of single harmonic ring loads to each of the nodes that may be in the 
tire-pavement contact region. 

TOTAL 3000 
RADIAL 
RING 
LOAD (lb) 

2000 


1000 


0 0.5 1.0 1.5 

RADIAL DEFLECTION (in) 

CROWN LOAD-DEFLECTION DATA CALCULATED WITH A UNIFORM RING LOAD APPLIED TO THE CROWN NODE 




SINGLE HARMONIC RING LOADS APPLIED TO A FINITE ELEMENT NODE 


Figure 2 
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TRANSFER- FUNCTION DEFINITION 

As a consequence of the linearity of the ring load-deflection response, the 
application of a single harmonic ring load produces a displacement field that varies 
circ umf erentially in the same harmonic as the applied ring load. Figure 3 gives the 
definition of the transfer function T^ as the ratio of the output and input 
amplitudes. Since each node responds differently, a transfer-function matrix, 

Tik] n> is needed to store the stiffness information generated by the ring loads. 


INPUT: Single Harmonic Ring Load 

OUTPUT: Single Harmonic Displacement 


cos n0 
cos n9 


TRANSFER FUNCTION ^ 


Tikin harmonic transfer function relating 

I displacement of node i to an n^*^ harmonic 
ring load on node k 


Figure 3 
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POINT LOAD VECTOR {p} AND THE DISCRETE FOURIER TRANSFORM (DFT) 

This application of the discrete Fourier transform uses an even number of 
points (N) , equally spaced around the circumference. The example shown in figure 
4 uses N = 8 points. A unit load is applied at any point, say point 0. The DFT 
of the load vector yields a set of N coefficients, G j , which are approximate values 
of the coefficients of the conventional Fourier series defined on the continuous 
interval 0 < 0 < 2 tt and representing the unit point load. The point load is 
applied, sequentially, in the radial, axial, and circumferential directions. 


INFLUENCE COEFFICIENT GENERATION 



I p } = { 1 . 0, 0, 0, 0, 0, 0, 0 J load vector 


1 it 

DFT G. = ^ E 9|cWj 
j N k=0 ® 


= e 

e 


•i2n/N 


9k = {Pf . Gj 


j = 0, 1. 


. N-1 


Figure 4 
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INVERSE DISCRETE FOURIER TRANSFORM (IDFT) AND THE INFLUENCE COEFFICIENTS 


Having the unit point load represented by a conventional Fourier series, whose 
coefficients a^ are approximately given by the DFT coefficients, the transfer func- 
tions are used, on each harmonic, to obtain the coefficients bn of the Fourier 

series representing the response of the nodal circle to the unit point load. The 
inverse discrete Fourier transform is then used to evaluate the displacements, Ujj^, 
at the N points. These displacements are the elements of the influence coefficient 
matrix as seen in figure 5 . 

INPUT SERIES COEFFICIENTS ^ 

OUTPUT SERIES COEFFICIENTS = ajik|n " ^ ^ik|n 


DFT OF DISPLACEMENT VECTOR 


IDFT S m = 0, 1, . . . , N-1 

i k 

INFLUENCE COEFFICIENTS j = 1, 2, . . ., 


shift: j = N 


SYMMETRY: = A.^,^ 


{d.j> - [A^j|(fhP|j£} 


Figure 5 



INFLUENCE COEFFICIENT MATRIX 


The Influence coefficient matrix relates the radial, axial, and circumferential 
components of the displacement of points on the tire surface to the radial, axial, 
and circumferential components of load at these points. The radial response parti~ 
tion shown in figure 6 is used to obtain a solution for frictionless contact, in 
which the axial and circumferential force components are known to be zero. The 
matrix here covers 3 points on each of 5 nodes. The point separation with this 
matrix is 11.25 degrees. 


= load at point i on node k 
d.j =• deflection of point j on node i 
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Figure 6 
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TOROIDAL SHELL CONTACT SCHEMATIC 


After the Inflation solution has been obtained, the tire model is deflected 
against a frictionless, flat surface. The contact surface is perpendicular to the 
wheel plane of symmetry and located at the specified loaded radius R^, as shown in 
figure 7. The vertical load and the contact boundary are unknown a priori. 



Figure 7 
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RADIAL DEFLECTIONS IN THE CONTACT REGION 


When the radius R£ is specified, the radial deflections are given approximately 
’^ij “ ^^i -1)A6] - R^, where R^ is the inflation radius of node i and 

A0 = 2tt/N is the point spacing. Since the contact half-angle is usually less than 
20°, the error in approximating the radial deflections by the above equation is not 
large. An initial estimate of the contact boundary is taken as the geometric 
intersection of the tire model and the contact surface. (See fig. 8.) 



Figure 8 
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LINE LOAD VECTORS 


The radial deflections within the contact boundary are known but the forces 
that produce these deflections are unknown. The tangential (axial and circumferen- 
tial) deflections within the contact boundary are unknown but the tangential forces 
are zero because the contact is frictionless. All surface forces are zero outside 
of the contact boundary. Since the number of unknowns (deflections and loads) is 
less than or equal to the number of equations established by the influence coeffi- 
cient matrix, an initial contact solution can be found. The contact boundary is 
then adjusted to exclude negative radial forces. Three to five boundary adjustments 
are normally needed to converge on the contact solution. Figure 9 shows the load 
vectors obtained in a solution for the G78-14 tire with 221 kPa (32 psi) inflation 
pressure. The elements of {p} are values of the line load at 32 points on the tire 
model equator. The other vectors give line load values in the right (and left) half 
of the contact region. Seven nodal circles are in the contact region in this ex- 
ample. 
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jpi} = |n.48, 5.67, 0, . . 0, 5.67 ( Ib/in 

{p2| = ( 16.27, 4.61 , 0, . . . , 0, 4.61 [ 

|p3| = 122.28, 0, . . ., 0( 

Ip**! = {13.78, 0, . . ., 0| 


Figure 9 











CONTACT PRESSURE DISTRIBUTION 


All of the contact pressure values (psi) calculated for the G78-14 tire with 
221 kPa (32psi) inflation pressure are shown in figure 11. The estimated location 
of the contact boundary is shown as a dashed oval. The contact boundary will be 
more accurately located when the density of points covered by the influence coeffi- 
cient matrix is Increased. The point density is limited only by the size and speed 
of the computer used to execute the tire model program. 
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CONTACT BOUNDARY 


Figure 11 
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AIRCRAFT TIRE SECTION 


The remainder of this presentation shows the contact solution calculated for 
a 32 X 8.8, Type VII aircraft tire. A theoretical and experimental study of this 
tire under inflation pressure loading was made by Brewer (cef • 5) . The photograph 
in figure 12 and the tire data shown in figure 13 are taken from reference 5. The 
white curve drawn on the tire section below marks the location of the carcass mid- 
surface. 
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MERIDIAN SECTION OF 32X8.8 AIRCRAFT TIRE 


Figure 12 
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MATERIAL PROPERTIES AND CARCASS GEOMETRY 


The parameters shown in figure 13 are used in a preprocessing subroutine to 
calculate homogeneous orthotropic properties for the finite element tire model. 


32 X 8.8 TYPE VII 
AIPXRAFT TIRE 

Material Properties and Carcass Geometry 


Rubber: Ep = 

450 psi, Vp 

= 0.49, Gp 

, = 151 psi 



Nylon Cord: 

= 156,000 

psi, = 

0.70, Gj, = 

700 psi 


Cord Diameter 

: d^ = 0.031 

in. 




Ply Thickness 

: h = 0.043 

in. (all plies) 



Cord Angle 6 

(measured from meridian) and Cord 

Density N, 

by Lift Formul 

Element 

B(deq) 

N(epi ) 

Element 

S(deq) 

N(epi ) 

1 

55.44 

25.42 

12 

46.97 

23.81 

2 

55.35 

25.39 

13 

45.34 

23.75 

3 

55.20 

25.34 

14 

43.71 

23.78 

4 

54.96 

25.26 

15 

41.99 

23.88 

5 

54.60 

25.15 

16 

40.25 

24.08 

6 

54.08 

25.00 

17 

38.64 

24.35 

7 

53.33 

24.79 

18 

37.13 

24.68 

8 

52.34 

24.55 

19 

35.73 

25.05 

9 

51.18 

24.32 

20 

34.39 

25.48 

10 

49.90 

24.10 

21 

33.37 

25.85 

11 

48.49 

23.93 




Construction: 6-ply, double bead 





Figure 13 
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FINITE ELEMENTS ON THE CARCASS MIDSURFACE 

The aircraft tire is modeled with 21 finite elements positioned along the car- 
cass midsurface. The ply structure in each element is specified separately in de- 
termining the homogeneous moduli for each element. Node 22 is a fixed node, posi- 
tioned to represent the tire bead which does not displace or rotate. (See fig. 14.) 


R 

(in) 



Figure 14 


LOAD VECTOR CORRECTION DATA 


The problem of calculating tire shape due to inflation pressure is highly 
nonlinear. As recognized by Stafford and Tabaddor (ref 6,), a successful solution 
can only be obtained by. a nonlinear finite element analysis which Includes updating 
the pressure load vector direction during the inflation solution procedure. Table 
1-1, in figure 15, gives the input load vector components, p^j and Pg, that are 
needed in order to have the resultant pressure load normal to the inflated tire 
model . 


CROWN 


DEFORMED MERIDIAN 


UNDEFORMED MERIDIAN 



TABLE 1-1. INPUT LOAD DATA FOR FINITE ELEMENT MODEL 
OF AIRCRAFT TIRE ANALYZED BY BREUER 


Element 

Number 

Rotation 

(deg) 

TIRELOAD Input Pressure (psl) 
P P 

1 

1.46 

94.97 

2.41 

2 

4.61 

94.69 

7.63 

3 

7.68 

94.15 

12.70 

4 

10.10 

93.53 

16.64 

5 

11.23 

93.17 

18.55 

6 

10.71 

93.35 

17.65 

7 

9.19 

93.78 

15.17 

8 

7.90 

94.09 

13.11 

9 

7.18 

94.26 

11.87 

10 

6.48 

94.39 

10.71 

11 

5.39 

94.58 

8.92 

12 

3.00 

94.87 

4.97 

13 

2.24 

94.93 

3.70 

14 

0 

95.00 

0 

21 

0 

95.00 

0 


Figure 15 
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CROVfN DISPLACEMENT VERSUS INFLATION PRESSURE 


The effect of correcting the load vector is clearly seen in figure 16. The 
finite element solution obtained when the pressure direction remains normal to the 
undeformed elements is indicated by A's. The solution found when the pressure is 
normal to the deformed elements is indicated by x's. This solution compares well 
with the calculation and measurements made by Brewer (ref. 5) . 



Figure 16 
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DEFLECTED MERIDIAN 


The deflected shape of the meridian passing through the center of contact is 
shown in figure 17 for the deflection 6 = 19 mm (0.75 in.). The tire load calcu- 
lated for this deflection is Fy = 9.76 kN (2194 lb). The distribution of contact 
pressure along the meridian is also shown. 



Figure 17 
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DEFLECTED EQUATOR 


The deflected shape of the equator and the distribution of contact pressure 
along it are shown in figure 18. Since only three points on the equator lie in the 
contact region, only a rough estimate of the circumferential location of the contact 
boundary can be made. 



Figure 18 
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TIRE LOAD VERSUS TIRE DEFLECTION 

Calculated values of tire load for specified tire deflections are shown in 
figure 19. 


(4 

/ 

/ /////// / 

32x8.8 TYPE VII / 

AIRCRAFT 

TIRE / 

95 psi INFLATION 

1 L 

1 1 


TIRE DEFLECTION 


Figure 19 






CONTACT AREA VERSUS TIRE DEFLECTION 

Calculated values of contact area for specified tire deflections are shown in 
figure 20. The calculated contact area increases in finite increments as additional 
points enter the contact region (as the load is Increased) . 



0 0.2 5 0.50 0.75 1.00 1.25 


TIRE DEFLECTION <T ( in.) 


Figure 20 


117 


CALCULATED CONTACT PRESSURE DISTRIBUTIONS 


All of the contact pressure values (psl) calculated for two different tire 
loads are shown in figure 21. The tire inflation pressure is 655 kPa (95 psi) . 


AIRCRAFT TIRE CONTACT PRESSURE DISTRIBUTIONS 



48 

✓ 

90 

21.3 

105 

34.4 

108 

28.3 

90 

34.4 

108, 

21.3 

105 


90, 

■s 



48, 



— ^ ^ 
^ 

,7 

V. 

S. 

,5 


,6 

21.3 

6 

34.4 

7 psi 

28.3 

6 

34.4 

6 

21.3 

5 


7 



(a) 6 = 0.75 in. , Fz = 2200 lb 



12.9 

130.4 

12.9 

V 

/ 




\ 

/ 

45.1 

187.9 

45.1 

\ 

/ 




\ 

/ 

57.3 

135.9 

57.3 

\ 

/ 




1 

\ 

42.9 

132.6 

42.9 

1 

1 

1 

38.2 

108.9 psi 

38.2 

I 

1 

1 

\ 

42.9 

132.6 

42.9 

1 

1 

\ 

57.3 

135.9 

57.3 

1 

\ 




1 

\ 

45.1 

187.9 

45.1 

/ 

\ 

12.9 

130.4 

12.9 

/ 

y 









" 



(b) 6 = 1 .00 in. , Fz = 3700 lb 


Figure 21 



EXTENSIONAL STRAINS DUE TO INFLATION PRESSURE 


The extensional (membrane) strains produced by inflation pressure only are 
shown in figure 22. In this plot, eg is the strain in the direction of the meridian 
and Bq is the strain in the circumferential direction. These strain distributions 
are axisymmetric and agree with the strains calculated by Brewer (ref. 5). 



0 1 2 3 4 5 6 7 


Crown Rim 

MERIDIAN DISTANCE S(in.) 


Figure 22 
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FORCE RESULTANTS DUE TO INFLATION PRESSURE 


The membrane forces (per unit length) produced by inflation pressure only are 
shown in figure 23. The forces Ng and Nq are in the meridional and circumferential 
directions, respectively. These force distributions are axisymmetric . 



Figure 23 
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ABSTRACT 


Finite element analysis is applied to study the large deflection and 
stress of a thin-walled pressurized torus loaded hy normal contact with a 
plane. The torus is found to have an elliptical footprint area, and con- 
siderable bulge occurs in the sidewall in the vicinity of the load plane. 
In large load ranges, the finite element calculations show compressive 
circumferential stress and negative curvature in the footprint region. 

An experimental study of the standing torus, using liquid metal strain 
gages , is outlined. Experimentally determined stresses are compared to 
those resulting from finite element analysis at various meridional and 
circumferential coordinates of the torus, including the footprint area. 
Circumferential strains compare favorably while meridional strains are 
higher in the finite element analysis, probably due to slippage of the 
boundary at the rim. 

This study utilized the STAGS finite element computer program. The 
purpose of the study was to evaluate the various program options for 
structural loading, for material modeling, for stress analysis, and for 
grid refinement. The experimental model, a thin-walled rubber tube mounted 
on a steel cylindrical rim, provided measured results to compare with 
various analytical trade-offs. It was found that there was almost no 
difference in predicted deflections or stress distributions between linear 
and nonlinear material description. However, the difference between lin- 
ear analysis and that of geometric nonlinearities utilizing incremental 
loading was marked . 


TORUS DIMENSIOMS AUD GLOBAL COORDINATE SET 


STAGS (ref. l) uses a first order shell theory to reduce a 3- 
dimensional structural prohlem to a dependency on 2 spatial coor- 
dinates , here the meridional coordinate ol and circumferential 
coordinate 0 . The deflection of the plane into the torus is 

5 . An example of an 11 x 15 torus grid is shown and dimen- 

sions of the cylindrical rim and inflated ruhher tube are given. 
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ORIGINAL PAGE IS 

OF POOR QUALITY 


EXPERIMENTAL SET-UP SHOWING ELLIPTICAL FOOTPRINT 

The tube was inflated to 0.03^ atm. (0.5 psi) and test loads of 
3h.3 N (7.7 lbs), 69.9 W ( 15.7 lbs), and 114.1+ N (25.7 lbs) were applied. 
Deflection, footprint area, circtmferential and meridional strain, and 
final pressure were observed. 




■ V. 




LIQUID METAL STRAIN GAGE AND CALIBRATION DEVICE 


In order to measure the large strain encoimtered , strain gages with 
a rubber housing and capillary column of mercury were used . The gages 
are inherently nonlinear by nature, and individual calibration is re- 
quired (ref. 2). Derivation of gage factor versus strain is given in the 




GAGE PLACEMENT AT 0° MERIDIONAL ANGLE - 
COMPLETE SYSTEM WITH DATA AQUISITION 



Calibrated gages were mounted on the torus with silicone rubber 
cement. The torus was inflated to test pressiire and rotated to the de- 
sired circumferential angle, load was applied, and readings were re- 
corded. The microprocessor shown below tabulated and recorded, on 
magnetic tape, data for strain versus circumferential angle at given 
values of load and meridional angle. 
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DATA ACQUISITION FLOW CHART 


The logic flow chart for the data acquisition scheme is shown helow. 
The buffer for the strain signals and the program controller with the 
stored logic of the acquisition system are encompassed in the micro- 
processor. The microprocessor is programmed in the BASIC language, and 
the digitizing of analog gage voltage is accomplished by a digital volt- 
meter. The three channels referred to are the meridional strain, the 
circumferential strain, and a dummy gage voltage. 
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LOAD-DEFLECTION CURVES FOR ANALYTICAL FINITE ELEMENT 

AND EXPERIMENTAL DATA FOR TORUS AT P. = 0.03^+ ATM 

1 

This figure shows points on the load-deflection curve obtained from 
a STAGS analysis using a 13 x 15 quarter torus at loads of 3^.3 N 
(T.7 Ihs), 69.9 N ( 15.7 ILs), and IIU.U N (25.7 Ihs ) . The deflections 
predicted from STAGS are lower since minimum potential theory is em- 
ployed and structural stiffness is overestimated with a coarse grid. 

As grid size is refined, predicted deflections approach the experi- 
mental values, but computer run time also increases (ref. 3). At a 
load of 3^.3 N, three different grid sizes were modeled. 
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■ 13x15 QUAETER TORUS GRIDS WITH 
. 3U.-3 W (T.T lUs) LOAD - NONLINEAR ANALYSIS 

This figure gives the computer-generated plot of the quarter torus with 
a 13 X 15 grid. Exploitation of the torus symmetry allows analysis of 
one-fourth of the structure, giving a grid density four times greater than 
that of a 13 X 15 grid applied to the full torus. The unloaded grid geom- 
etry is overplotted . on the the deformed geometry. Note the plane surface 
at the load point (footprint) and the side wall hnilge. The nonlinear 
geometric analysis option of STAGS gave convergence to a final deflection 
of 6k% of the experimental value while the linear analysis yielded only 
21^ of the experimental value. This result is significant when compared 
to toroidal analysis using other finite element procedures (refs. 4 and 5). 




CALCULATED STRESS DISTRIBUTION AT 90 ° MERIDIONAL 
ANGLE FOR 13 x 15 QUARTER TORUS GRID - L = 3^^.3 N (7-7 lbs.) 

Calculated values of meridional stress and circumferential stress 
versus circumferential angle 0 are shown below. The variation of 
stress in the footprint and bulge regions is significant. The values 
remain constant away from the contact area and are mainly due to in- 
ternal pressure. These values, with this grid refinement, agree 
favorably with other investigators (ref. 6). 





CALCULATED MERIDIONAL STRESS DISTRIBUTION, , 13 x 15 
QUARTER TORUS GRID - NONLINEAR ANALYSIS, LOAD = 3h.3 N (T.7 ILs) 

The three dimensional plot below shows meridional stress dis- 
tribution over the surface of the torus at a 3^.3 N load. Even 
tensile stress appears in both the meridional and circumferential 
directions in the upper half of the torus , as a result of internal 
gas pressure. However, the compressive effect of contact with the 
flat plate reduces these tensile stresses in the footprint region. 
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CALCULATED CIRCUMFERENTIAL STRESS DISTRIBUTION, , 13x15, 
QUARTER TORUS GRID - NONLINEAR ANALYSIS, LOAD = 3^.3 N (j.j lbs) 

A three dimensional plot of circumferential stress distribution 
over the upper half of the torus is shown below. As with the previous 
figure, uniform tensile stress due to gas pressure is significantly 
reduced in the footprint region, to almost 0 N/mm^ in the center. Results 
of an intermediate load step of 69-9 N (15.7 Ihs) are contained in 
the literature (ref. T). 
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13 X 15 QUARTER TOBUS GRID WITH 
1lU.i+ N (25.7 lbs) LOAD - NOHLIHEAE ANALYSIS 

The figure below shows the STAGS generated plot of a quarter 
torus in the unloaded and deformed configurations under a test 
load of 11i^.i+ N. The deflection has increased significantly and 
the degree of bulge in the sidewall is more pronounced. The center 
of the footprint region shows negative curvature. 


CALCULATED CIRCUMFERENTIAL STRESS 
DISTRIBUTION, , AT LOAD = IlUA N (25. T ILs . ) 

A three dimensional plot of circvunferential distribution over 
the quarter torus is given below for the IIU.I+ N load. Note the 
scale change from that of the previous two three dimensional plots . 
Increased compressive stress in the footprint region is seen to cause 
a stress reversal in the center of the footprint region. The stress 
becomes compressive near the center, then suddenly reverses and 
becomes tensile, suggesting that a local limit point on the structure 
load deflection path has been reached and snap-through has occurred. 


3.5 X 10"2 U/m^ (50 PSD 



fOOTPRIHT REGION 



FOOTPRINT OF AN INNER TUBE SHOWING NEGATIVE CURVATURE 




The figure below is a photograph of an inner tube under rela- 
tively heavy load. An initial eccentricity has been introduced. 
The negative curvatirre near the center of the footprint is consid- 
erably reduced in the presence of friction at the contacting sur- 
faces (not modeled here). 
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Strain (percent) 


EXPERIMENTAL AND CALCULATED STRAINS AT 
-10° MERIDIONAL ANGLE - LOAD = 3^.3 N (j.j lbs.) 

The figiire below- compares measured and calculated meridional 
and circumferential strains for varying 0 at « = -10°. This -10° 
meridional angle places the strain gages in the region near the rim. 
Good agreement is achieved in obsemred and calculated values of 
circumferential strain. The disagreement between observed and cal- 
culated meridional strain is attributed to slipping of the membrane 
iinder the rim on the test fixture . The computer model had a fixed 
boundary at the rim. 


B Meridional Strain (Predicted) 
Meridional Strain (Measured) 

^ Circumferential Strain (Predicted) 
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EXPERIMENTAL AND CALCULATED STRAINS AT 
90° MERIDIONAL ANGLE - LOAD = 69-9 N (15.T lbs.) 

This figure compares measured and calculated strains for vary- 
ing circumferential angle 6 at a = 90°. At this farthest point 
from the rim, meridional strains are in good agreement. Note here 
that for all measurements below a = 30°, the strain gage is in the 
footprint region. Measiired strain did not become compressive here, 
probably due to friction between the contacting surfaces. Reports of 
strain comparison at other loads are contained in the literature 
(ref. 8) . 
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APPENDIX 


In describing the operation of a liquid metal strain gage, one 
assumes that the liquid metal column obeys the fundamental resist- 
ance relation (there are no voids in the column). Further, the 
resistivity of the conducting medium is assumed constant, and the 
volume of the capillary cylindrical cavity is assumed constant. 
Rubber is essentially incompressible; the assumption of constant 
volume under elongation is a valid one. The fundamental resistance 
is 


R _ CL _ p^ . 
A V 


( 1 ) 


when p = resistivity, ohm-in. 

L = length, in. 

A = cross-sectional area, in. 
V = volume, in. = aL 


Using a Taylor series expansion about the initial length, 
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( 2 ) 


where ^ - strain, in. /in. 

Equation 2 establishes a quadratic relationship between resist- 
ance change and strain. It is seen from ( 2 ) that for small strains, 
the gage factor is two. 
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EXPERIMENTAL AND FINITE ELEMENT STUDY OF A 
STANDING TORUS UNDER NORMAL AND TANGENTIAL LOADS 


Donald R. Flugrad and Bruce A. Miller 

Iowa State University 


EXPANDED ABSTRACT 

There continues to exist a strong interest in the load deflection attributes of 
both automotive and aircraft tires. Considerations of safety, handling and perform- 
ance during such maneuvers as landing, braking and cornering are Intimately tied to 
the adequacy of a vehicle's tire design. In the past, the qualitative and quantita- 
tive information needed for a basic understanding of tire characteristics has been 
derived by strictly empirical methods. Of late, however, numerical techniques such 
as finite element analysis have been brought to bear on certain aspects of the 
problem. 

This paper describes work that is presently in progress at Iowa State University 
to determine the effect of a combined external load consisting of a normal component 
and a tangential braking force applied to an inflated torus. Experimental results 
obtained by photographic study of the contact area between the torus and load plate 
are presented as well as measurements of the vertical and horizontal displacement of 
the torus under load, A numerical procedure for displacement analysis is developed 
in which the finite element program STAGS (Structural Analysis of General Shells) is 
used in an Iterative manner to produce a flat, horizontal footprint surface under 
force loading. The redistributed force distribution obtained by the iterative 
process is displayed along with computed meridional and circumferential stresses. 
Finally, an extension of the iterative method is introduced which eliminates the need 
to experimentally determine the footprint area. 
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LINE DRAWING OF TEST FIXTURE 


Figure 1 shows the torus mounted on a rim and axle. The axle is clamped to the 
frame to prevent rotation of the unit. Two load tables work independently to produce 
an externally applied load on the torus. The table directly below the torus acts 
through the pulleys at the top of the fixture to pull the plate in contact with the 
torus upward, providing a normal load. The table shown to the left of the fixture is 
then loaded to exert a tangential force in the fore-aft of braking direction. 






Figure 1 







PHOTOGRAPH OF TEST FIXTURE 


ORIGINAL PAGE IS 
OF POOR QUALITY 


The actual test fixture is depicted in Figure 2 with weights mounted on the 
normal load table directly below the 8.00 x 20 inner tube used as the torus. A 
plastic milk container is shown in place of the tangential load table of Figure 1. 
By introducing water into the milk container, it is possible to vary the tangential 
load continuously instead of incrementally as would be the case if weights were 
placed on a load table. The plate in contact with the torus is made of transparent 
Incite allowing photographs of the footprint area to be taken from below the plate. 
Graph paper "targets" are affixed to the Incite plate for the purpose of tracking 
the vertical and horizontal displacement of the plate under prescribed load condi- 
tions. The manometer at the right is used to monitor the internal pressure of the 
torus. 



Figure 2 


143 


ORIGINAL PAGE IS 

PHOTOGRAPHS OF FOOTPRINT AREAS POOR QUALITY 


The two photographs of Figure 3 show footprint areas viewed through the trans- 
parent Incite plate. In both cases the torus is subjected to a normal load of 70 N 
and an internal pressure of 3450 Pa. Figure 3(a) was taken with no tangential load 
applied, and figure 3(b) shows the footprint area produced when a tangential 
load equal to 80% of the normal load is added. The photographs show that the contact 
area for each case is generally elliptical in shape. Closer inspection of the two 
photos verifies that there is little if any distortion introduced by the inclusion of 
a tangential component of force acting on the torus. 









(b) Load = 70 N; braking 








EXPERIMENTAL MEASUREMENT OF DISPLACEMENT 


The experimental set-up for the measurement of vertical and horizontal displace- 
ments of the torus under prescribed load conditions is depicted in Figure 4. In 
actuality the center of the inner tube on its rim and axle is held stationary and the 
displacement of the Incite plate in contact with the torus is measured. Three graph 
paper "targets" are attached to the load plate, two of which are visible hanging from 
the front edge of the plate in the photograph. The third is attached to the back 
edge of the plate. Three separate optical devices are trained on the grid patterns 
of the "targets" to record vertical and horizontal displacements of three distinct 
points on the plate as normal and tangential braking forces are applied to the torus. 
Only two of the three optical devices appear in the photo. The third, a surveying 
transit with a much longer focal length, is stationed some distance from the test 
fixture. The three separate measurements of displacement are then used to estimate 
the vertical and horizontal movement of a point located at the center of the foot- 
print area. 






VERTICAL DISPLACEMENT OF TORUS 


The graph of Figure 5 illustrates the relationship between the vertical dis- 
placement of the torus and the magnitude of the applied tangential braking force for 
a constant normal load of 70 N. Experimental results plotted on the graph are those 
obtained by the set-up of Figure 4- The solid line represents a least squares linear 
fit of the experimental data. Also shown are results obtained from STAGS, a finite 
element program particularly suited for analysis of plates and shells, for three 
different tangential load cases. The tangential force associated with the horizontal 
axis is a non-dimensional quantity which represents the braking force as a decimal 
fraction of the applied noirmal load. It is evident that the vertical displacement of 
the torus is relatively insensitive to the tangential force applied. It is also clear 
that the finite element analysis predicts a vertical displacement which is consist- 
ently less than that obtained by measurement. 


VERTICAL DISPLACEMENT. 70 N NORMAL LOAD 
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HORIZONTAL DISPLACEMENT OF TORUS 


The graph of Figure 6 is similar to that of Figure 5 except that horizontal 
displacements of the torus are recorded rather than vertical displacements. A least 
squares cubic fit curve was chosen to best represent the experimental data after 
first considering both linear and quadratic forms. Of the three cases analyzed by 
finite element techniques, the first two are seen to agree remarkably well with the 
experimental results. For the greatest tangential force, however, STAGS predicts a 
horizontal displacement which is approximately only 50% of the measured value. Such 
a large discrepancy suggests that some undetected slip may have occurred between the 
lucite plate and inner tube for this case involving a high tangential braking force. 
Separate consideration of the three results arrived at by finite element analysis 
indicates the possibility of a more nearly linear relationship between horizontal 
displacement and tangential force. However, more data would have to be gathered 
before such a conclusion could be drawn. 


HORIZONTAL DISPLACEMENT. 70 N NORMAL LOAD 



Figure 6 
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FINITE ELEMENT COORDINATES 


Nominal dimensions of the inner tube used in this study are shown in Figure 7 
along with the two coordinate angles required by STAGS to describe a toroidal shaped 
element. The angle 6 is the circumferential angle. It varies from -180° to +180° 
for the torus. The meridional angle, a, varies from -42° to 90°. At a = -42° the 
torus is assumed to be clamped at the rim. Because there exists a plane of symmetry 
through the equator of the torus with respect to geometry as well as to the externally 
applied loads, it is possible to perform a complete finite element analysis by con- 
sidering just half the inner tube. Hence, the maximum value for a is just 90° rather 
than 222°. 


DIMENSIONS AND COORDINATES 



SECTION A-A 


Figure 7 
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FINITE ELEMENT DISPLACEMENT RESULTS 
FOR QUARTER TORUS WITH NORMAL LOAD 


Shown in Figure 8 is the deformed as well as the undeformed geometry of the 
quarter torus analyzed by Mack (ref. 1) in a previous study. The finite element 

analysis was conducted with STAGS using a 13 x 15 grid pattern to define the loca- 

tion of the element nodes. A quarter torus was adequate in this case because there 
are two planes of symmetry in the absence of tangential loading. The 70 N normal 
load was evenly distributed over the nodes located within the experimentally 
determined footprint area. As can be seen, the evenly loaded nodes situated along 
the bottom of the torus in both views appear to lie on a relatively flat horizontal 

plane. Since the load plate is assumed to be horizontal and rigid, this result 

justified the use of an even load distribution in the footprint area for this 
particular case. 
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INITIAL RESULTS FROM FINITE ELEMENT 
DISPLACEMENT ANALYSIS WITH TANGENTIAL LOADING 


The deformation of the half torus used to analyze the effect of combined normal 
and braking forces is illustrated in Figure 9. A 10 x 21 grid pattern was chosen in 
an attempt to provide a fine enough mesh for reasonably accurate results at a 
reasonable cost. The 70 N normal load is identical to that used for the quarter 
torus of Figure 8. In this initial run both the normal and tangential forces were 
evenly distributed over the nodes experimentally determined to lie within the foot- 
print area. Unlike the normal load case of Figure 8, however, the loaded nodes along 
the bottom of the torus do not appear to lie on a horizontal plane. The tangential 
force, which acts to the left in the left-hand view, has the effect of tilting that 
portion of the toroidal surface meant to be in contact with the horizontal load plate. 
In an effort to correct this problem, an iterative method was devised to achieve a 
comparatively flat and horizontal surface across the bottom of the torus by systemati- 
cally redistributing the forces on the nodes within the footprint area. 


HALF TORUS, 70 N NORMAL, 801 TANGENTIAL, INITIAL RUN 



Figure 9 
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REDISTRIBUTION OF FORCES 


The simplified model of Figure 10 illustrates the method developed to redistrib- 
ute the forces in the footprint region. Suppose several springs of different lengths 
with different spring rates are suspended from a horizontal beam as shown. The total 
force to be exerted upward on the group of springs is known, and the individual 
forces to be applied to each of the springs is to be determined so that they all 
deflect to a common vertical equilibrium position denoted by Zg. An arbitrary dis- 
tribution of forces is chosen to begin the procedure. For example, the total force 
might initially be uniformly distributed over all the springs. The resulting 
deflection of each of the springs is observed, and a corresponding set of spring 
rates is calculated. From a consideration of the total force, the common equilibrium 
position, Zg, is determined from the equation given in the lower right hand corner of 
the figure. Once Zg is known the forces to be exerted on each of the individual 
springs are calculated from the equation provided for Fj. For a linear set of springs 
which act independently, these forces will, in fact, produce a common equilibrium 
position for all the springs. The situation with the inflated torus, however, is 
more complicated and requires an iterative scheme. The nodes within the footprint 
region are first loaded evenly. The STAGS analysis then gives the resulting deflec- 
tions from which individual nodal spring rates are determined. A trial value is next 
computed for the equilibrium position just as in the case of the simplified model, 
and forces to be applied to the individual nodes are evaluated. These forces are 
then used to perform another STAGS analysis. Since the spring rates of the nodes are 
expected to be nonlinear, and since the nodes do not act independently, this second 
application of STAGS will not, in general, produce a flat, horizontal surface on the 
torus. Therefore, the entire procedure must be repeated until the deflected nodes 
are found to lie on a surface which is sufficiently flat and sufficiently horizontal. 



F, = K,(Z, -Z,) 






Figure 10 
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FINITE ELEMENT DISPLACEMENT ANALYSIS 
FOR FIRST ITERATION 


After redistributing the forces on the torus of Figure 9 by the method described 
in Figure 10, the STAGS analysis was performed again. The results of that analysis 
are shown in Figure 11. As can be seen, the footprint area is still not flat, but 
the severe tilt evident in Figure 9 has been greatly reduced. 


HALF TORUS, 70 N NORMAL, 80% TANGENTIAL, FIRST ITERATION 



Figure 11 



FINITE ELEMENT DISPLACEMENT ANALYSIS 
FOR FINAL ITERATION 


Figure 12 depicts the displacement results for the third iteration on the torus 
of Figures 9 and 11. Although some waviness is still apparent in the footprint 
region, the surface is relatively flat and horizontal. In fact, further iterations 
failed to produce any significant improvement. The common equilibrium position 
calculated for this case was used to establish the vertical deflection plotted in 
Figure 5 for a non-dimensional tangential force of 0.8. Only the normal components 
of force were considered in the redistribution process. Once the normal forces were 
established for each of the nodes, a tangential component of force equal to 0.8 of 
the corresponding normal force was added. Consequently, the resulting horizontal 
displacements of the nodes in the footprint region showed some variation. In order 
to arrive at a single value to characterize the horizontal deflection of the torus, 
the average of the horizontal displacements of the nodes across the center of the 
footprint area at a circumferential angle of 0° was computed. This value was then 
plotted in Figure 6 as one of the three data points found by the finite element 
method. 

HALF TORUS. 70 N NORMAL. 80Z TANGENTIAL. FINAL ITERATION 



Figure 12 
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FORCE DISTRIBUTION FOR FINAL ITERATION 


The footprint force distribution determined for the final Iteration of Figure 12 
is displayed in Figure 13. The sense of the tangential force applied to the torus 
is down and to the left with respect to the three dimensional plot. As shown, the 
force is greater along the sides of the footprint and at the left hand leading edge 
than in the middle. While the distribution , appears to vary significantly, the 
minimum force at the middle of the footprint is actually 83% of the maximum force 
at the leading edge. 


FOOTPRINT FORCES/ 70 N NORMAL 80% TANGENTIAL FINAL ITERATION 



Figure 13 
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MERIDIONAL STRESSES FOR INITIAL RUN 


The meridional stresses computed by STAGS for the Initial run of Figure 9 are 
shown in Figure 14. The footprint area extends from approximately -25° to +25° in 
the circumferential direction and from approximately 90° to 45° in the meridional 
direction. The three dimensional display suggests the existence of a great deal of 
variation in the stresses throughout the torus with the stress even becoming compres- 
sive for a small area in the neighborhood of 0 = 0° . However, caution must be 
exercised in using this plot and those of the next three figures. There have been no 
experimental stress measurements made to verify these results, and no other grid 
patterns have as yet been tried to determine the effect of changing the finite ele- 
ment mesh on the computation of stresses. Another point that should be made with 
respect to Figures 14 through 17 is the fact that in executing the STAGS analysis the 
nodes located at circumferential angles of -180° and +180° were constrained to remain 
in a vertical plane. Hence, the stresses calculated at these extreme values of 0 
turned out to be quite different even though they represent the same points on the 
torus. Elimination of this constraint could have a significant effect on the pre- 
dicted stress distribution. 


MERIDIONAL STRESS. 70 N NORMAL. 80% TANGENTIAL. INITIAL RUN 



Figure 14 
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MERIDIONAL STRESSES FOR FINAL ITERATION 


The meridional stresses for the final iteration of Figure 12 are depicted in 
Figure 15. Comparison with Figure 14 indicates that redistribution of the forces in 
the footprint area produces a dramatic effect on the stress distribution. 


MERIDIONAL STRESS. '70 N NORMAL/ 80% TANGENTIAL. FINAL ITERATION 



Figure 15 
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CIRCUMFERENTIAL STRESSES FOR INITIAL RUN 


Figure 16 shows the circumferential stresses obtained from the STAGS analysis 
for the initial run of Figure 9. It is interesting to note the sharp decrease in 
stress at the leading edge of the footprint at a circumferential angle of approxi- 
mately -40° . This seems reasonable in light of the fact that the tangential force 
is directed to the left with respect to the three dimensional plot. 


CIRCUMFERENTIAL STRESS. 70 N NORMAL. 80% TANGENTIAL. INITIAL RUN 



Figure 16 
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CIRCUMFERENTIAL STRESSES FOR FINAL ITERATION 


Figure 17 depicts the circumferential stress results for the final iteration of 
Figure 12. The high stresses in the footprint region and the sharp decrease just to 
the left of the contact area are even more pronounced for this case with redistrib- 
uted forces than for the initial run of Figure 16. 


CIRCUMFERENTIAL STRESS. 70 N NORMAL. 80% TANGENTIAL. FINAL ITERATION 



Figure 17 
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FINITE ELEMENT DISPLACEMENT RESULTS 
FOR INITIAL RUN OF SIX NODE CASE 


In all the finite element work described thus far, the area of the footprint 
region was known beforehand from photographs taken on the test fixture of Figure 2. 
Nodes that were found to fall within the footprint area were then subjected to normal 
and tangential forces. In an attempt to eliminate the need to experimentally 
determine the contact area, a procedure was developed to start the iteration process 
described in Figure 10 by initially loading just six nodes centered about 0=0° and 
a = 90°. Theoretically, it should be possible to start with a single node, but it 
was decided that such a choice would only delay convergence to a final solution. The 
displacement results obtained when the full normal load of 70 N and the 80% tangen- 
tial braking force were uniformly distributed over just six nodes are displayed in 
Figure 18. After spring rates were calculated for the six nodes and the first esti- 
mate of the common equilibrium position was determined, it was discovered that 102 
nodes lay below Zg. In the next run the total normal and tangential forces were 
evenly distributed over all 102 nodes. The next estimate of Zg turned out to be 
considerably lower with only 83 nodes located below it. A second iteration reduced 
to 60 the number of nodes to be loaded. This corresponds to the number of nodes 
found experimentally to be within the footprint area for this particular load case. 
Three additional iterations were conducted to redistribute the forces by the method 
of Figure 10 in order to refine the flatness of the contact area. 

SIX NODE CASE. HALF TORUS, 70 H NORMAL, SOX TANGENTIAL, INITIAL RUN 



Figure 18 
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FINITE ELEMENT DISPLACEMENT RESULTS 
FOR FINAL RUN OF SIX NODE CASE 


The deformed geometry for the fifth iteration of the six node case is shown in 
Figure 19. Although not perfectly flat, the results are comparable to those of the 
final iteration depicted in Figure 12. While not included in this report, the final 
force distribution was also found to be very similar to that shown in Figure 13. 


SIX NODE CASE. HALF TORUS. 70 N NORMAL. 80Z TANGENTIAL. FINAL RUN 



Figure 19 
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OVERVIEW OF NASA TIRE EXPERIMENTAL PROGRAMS 


John A. Tanner 

NASA Langley Research Center 
ABSTRACT 

The experimental measurement of various tire properties is an important step in 
the tire modeling process. Some properties are needed to define the complex loading 
distributions that must be imposed on a tire model to simulate typical operational 
environments. Others, such as tire stresses and displacements resulting from 
various loading conditions, can he used to assess the accuracy of candidate finite 
element codes. In many cases these properties are used to establish the tire 
response characteristics which are useful in the simulation of operational problems 
and support the design of vehicle suspension systems and aircraft landing gear sys- 
tems . 


This paper reports on the interim results from a number of ongoing aircraft 
tire experimental programs. These programs are designed to measure profile growth 
due to inflation pressure and vertical loading, contact pressures in the tire 
footprint, and a number of tire mechanical properties Including spring, damping, 
and relaxation characteristics. 



APPARATUS FOR TIRE PROFILE MEASUREMENTS 


of Sd 'S 

POOR QUALITY 


Profile growth measurements of a 40 x 14 aircraft tire for various inflation 
pressure and vertical loading combinations are under way using the equipment shown in 
figure 1. The tire is mounted in the tire vibration stand described in reference 1 
and a pair of linear variable displacement transducers is connected to an x-y 
plotter to produce full scale profiles. Figure 2 is a close-up view of the tire 
showing the markings that were used to insure repeatable test results. 



Figure 1 



Figure 2 
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PROFILE GROWTH DUE TO INFLATION PRESSURE 


Typical profile growth data for the 40 x 14 aircraft tire during inflation 
are presented in figure 3. For this partictilar test the tire was inflated to 
1.17 MPa (170 psi) from an uninflated state and the arrows in the figure denote 
the direction of growth for selected points on the tire periphery. The profile 
growth involved both radial and tangential tire displacements. The tire inflation 
problem is axisymmetric and the tire displacements are symmetric about the 
tire equator. 
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PROFILE GROWTH DUE TO COMBINED INFLATION PRESSURE AND VERTICAL 

LOADING 

Typical growth measurements of a 40 x 14 aircraft tire for vertical loads 
ranging up to 89 kN (20 000 Ibf) are presented in figure 4. For this test the 
inflation pressure was 1.07 MPa (155 psl) and the profile measurements were along 
a line which was displaced 5° from the contact centerline as shown in the sketch. 
The ground lines for each of the vertical loading conditions are also shown in 
the figure. The arrows again denote the direction of growth and in this case 
depict the growth of the sidewall bulge for increasing vertical load. The dis- 
placement data from these loading conditions should exhibit symmetry about the tire 
equator but the axial symmetry associated with the Inflation process is now lost. 
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TIRE FOOTPRINT 


ORfGfNAL PAGE (S 
OF POOR QUALITY 


The forces and moments developed in the tire footprint must be modeled 
accurately to simulate the tire operational environment. Figure 5 xs a tyPi^al 
footprint from a 40 x 14 aircraft tire inflated to 1.07 MPa (155 psx) and subjected 
to a vertical load of 89 kN (20 000 Ibf). Instrumentation is currently under 
development to measure the pressure distribution and tangential friction force 
variation throughout the tire/surface contact region. 



Figure 5 
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FOOTPRINT LENGTH 


The length of the tire footprint for two tire sizes is presented as a 
function of tire vertical deflection 6 in figure 6 where both parameters have been 
nondimensionalized by the tire outside diameter d. The data for the two tire sizes 
can be faired in a least squares manner by a single-valued, nonlinear curve defined 
by the following expression (ref. 1), 


Lj/d = 1.66 


(6/d) - (6/d)^ 


( 1 ) 


This expression is similar to an expression in reference 2 for various type 
VII aircraft tires. Also included in the figure is the curve for the expression 
which defines the relationship between footprint length and vertical deflection if 
the tire is not distorted by the vertical load: 


L^/d = 2 


'(6/d) - (6/d)' 


( 2 ) 


Without distortion the footprint length equals the length of the geometric 
chord formed by the intersection with the ground plane of a circle having a 
diameter equal to that of the tire. The data for both tires indicate that the 
footprint is about 83 percent of the geometric-chord length. 



Figure 6 
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STATIC LOAD-DEFLECTION CHARACTERISTICS 


Typical lateral load deflection curves for an aircraft tire are presented in 
figure 7. This family of curves was generated by varying the lateral load while 
holding the vertical load constant. For this test the lateral loading was decreased 
with each successive cycle as identified by the numbers in the figure which denote 
the loading sequence. The tire response to this loading condition is characterized 
by a number of large hysteresis loops. The area enclosed within each loop is a 
measure of dissipated energy or hysteresis loss and the slope of the load deflection 
curve is a measure of the tire spring characteristics. Similar load deflection 
curves can be generated in the vertical and fore and aft or braking directions to 
establish additional tire response characteristics. 



Figure 7 
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HYSTERESIS LOSS 


The static lateral hysteresis loss of a tire can be obtained from load 
deflection curves such as that shown in figure 7 by measuring the area enclosed 
within the lateral loading and unloading cycle. Typical values of the lateral 
hysteresis loss ratio are plotted as a function of lateral load in figure 8. The 
lateral hysteresis loss ratio, the ratio of hysteresis loss to total energy input 
to the tire, is obtained by dividing the area within the hysteresis loop by the 
area under the load deflection curve including the hysteresis loop. The data in 
the figure indicate that increasing the lateral load Increases the energy input 
lost to hysteresis, and the influence of vertical load is inconsequential. 



Figure 8 
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STATIC AND DYNAMIC LATERAL SPRING RATES 


Typical tire lateral spring rates obtained from static and dynamic tests are 
plotted in figure 9 as a function of vertical load. As shown in the sketch, two 
spring rate values were obtained from the static tests. One spring rate corresponds 
to the linear portion of the load-deflection curve during load application and is 
denoted by the short dashed line in the figure. The other static spring rate 
approximates the slope of the load-deflection curve at Initial load relief and is 
denoted by the long dashed line. These two static spring rates represent the 
lower and upper bounds for each vertical loading condition. The dynamic spring rates 
were determined from free vibration tests described in reference 3 and are denoted 
by the solid line. The data presented in the figure indicate that dynamic spring 
rates lie between the two bounds established by the static tests. 


SIZE 49 X 17 TIRE 
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LATERAL DEFORMATION OF TIRE PERIPHERY 

Another parameter available from static tests is the static relaxation length 
which is obtained from displacements of the tire equator at various angular positions 
around the circumference. Such measurements taken on a 49 x 17 aircraft tire under 
four combinations of vertical loading Fg and lateral loading Fy are presented in 
figure 10. Lateral deformation of the tire equator \ is plotted as a function of 
the tire peripheral angle 6 measured from the footprint centerline. The lateral 
displacement of the free-tread periphery of the tire near the edge of the footprint 
varies exponentially with circumferential angular position and can be expressed 
in the following form: 

X = Ae “®^^s (3) 

where A is a constant, s is the circumferential distance from the footprint edge, 

and L is called the tire static relaxation length, 
s 
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Figure 10 




VARIATION OF STATIC RELAXATION LENGTH WITH 
VERTICAL DEFLECTION 

Relaxation lengths computed from data such as that presented In figure 10 are 
plotted in figure 11 as a function of tire vertical deflection. Both parameters 
are nondimensionallzed by the tire width w. Data from two tire sizes are presented 
and faired by separate linear curves. Relaxation lengths from both tires tend to 
decrease with increasing tire deflection. 
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CONCLUDING REMARKS 


Some of the data reviewed in this paper are typical of the experimental measure- 
ments necessary to verify the performance of tire analytical models and to establish 
footprint loading conditions for tire contact problems. Other data presented are 
used in defining tire spring and damping properties and the transient response 
characteristics necessary to model the tire's role as a component of a vehicular 
suspension system or an aircraft landing gear system. Experimental measurements of 
tire properties will continue to be an Important step in the tire modeling process, 
regardless of the level of sophistication of the modeling technique. 
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ADHESION TESTING OF AIRCRAFT TIRES 


Stephen N. Bobo 

U.S. Department of Transportation 
Transportation System Center 


In December 1979 the FAA issued a new Technical Standard Order TSO-C62c to 
all users and manufacturers of aircraft tires. It was designed to upgrade the 
testing required to meet minimum airworthiness standards. 

These changes to the testing requirements for new tires necessitated similar 
Improved standards for retreads used in the national air carrier fleet. 

Accordingly, an advisory circular (ref. 1) was prepared for comment which 
upgraded the testing standards for retreads to reflect the changes made in testing 
new tires under TSO-C62c. The advisory circular recommending the new dynamometer 
testing requirements called for testing every retread level of every tire size 
in an effort to accumulate sufficient tests and data to provide confidence that 
the retreading process including casing selection contained procedures which would 
provide for the continued airworthiness of the tires in service. 

However, the number of tires to be tested to accumulate confidence would have 
presented an unacceptable and unrealistic cost to retreaders and their customers 
and an alternative approach was necessary. 

For many years tire manufacturers and retreaders have been using laboratory 
adhesion tests as means for determining the effectiveness of the vulcanizing pro- 
cess in adhering the various tire components to one another. Adhesion testing 
appeared to offer a less burdensome alternative to replace some of the dynamometer 
tests recommended in the AC, Accordingly, test results and data were requested 
from retreaders who had used adhesion testing. 

All of the American retreaders of aircraft tires submitted data, as did 
Goodyear, which obtained additional data from European adhesion tests. For some tires 
the Navy has required adhesion tests as a part of their purchasing procedure and 
this data was also made available. Such data was collected from various sources 
for over 700 tires, both commercial and military. 

In meetings with Industry, the FAA was presented with the industry consensus 
regarding the use of ply and tread adhesion tests to qualify tire retreading 
process specifications. The FAA has accepted this means of testing as one which 
can be used in the qualification of a tire retreading process specification. 

The adhesion testing procedure used by most retreaders was a modification of 
the Hascar-Reiger method (ref. 2) in which a 1" strip of rubber is slit and intro- 
duced into a tensilometer . Variations of this method are described in references 
1 and 3 to 6. 

Figure 1 gives an example of the output from an adhesion test. The data is 
subject to wide variations in Interpretation because of the stress-tear-relaxation 
characteristic of the rubber sample. Several methods of reporting the data have 
been adopted, including averaging all maximum values, averaging all minimum values 
and taking the mean of the maximum and minimum average values. The reporting 
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sources usually list the method of recording and this is so noted in the data. 
Several other variations in method account for variations in test results, the 
three most important of which relate to sample preparation. Some laboratories 
attempt to cut the sample to the exact dimension, others correct for errors in 
size by normalizing the cut dimension to 25*4 mm (1 in.) width. Some laboratories 
cut the sample approximately 3.2 mm (1/8 in.) oversize and using a razor blade 
slit the intended path of travel of the tear line around the edge of the sample. 

Several Individuals have reported high values of adhesion when excess sample 
rubber thickness is not cut away, however laboratories as a rule do nothing to 
alter the thickness of the sample. 

The location of the tear region varies depending upon the agency requesting 
the testing, and when it is known this information has been included in the data. 
Most organizations have reported buff line adhesion data, although some have 
reported maximum rather than average results. 

Some organizations have reported adhesion data from the outside of the outer 
ply while some have reported between the second and third ply. Since there was 
no statistically detectable difference between these reported values they were 
lumped in the data. 

The data was tabulated and placed in a data base called BANK (ref. 7). The 
fields are described in a listing. Most of the data on buffline adhesion is taken 
from TAV although TMX contains some buffline data. These two fields were 
separated because of uncertainty about the method of reading the primary recorder 
traces . 

PAV gives values of outer ply adhesion. Some readings of maximum and minimum 
adhesion averages were available and these were recorded as PMX and PMN. Tire 
size is structured so that mathematical transformations such as linear regression 
or rank order correlation can be performed on the size variable. 

Other information, such as R level, durometer tensile, and elongation measure- 
ments, is included where a sufficient amount of data was obtained. 

The BANK program provides an interesting first level statistic printout of 
the data in each field (Figures 2 through 13). These include mean, standard 
deviation, and maximum and minimum values, as well as a data histogram. 

In order to use simplified procedures for establishing minlmiom adhesion 
thresholds and realistic test sample sizes it is Important to confirm the character 
of the distribution of the data. This may be accomplished by analysis of the data 
in terms of the probability that it fits on a normal distribution curve (Figure 
14). 


Figure 15 is a scatter plot of the probability that any given sample will lie 
below a given value, against ordinal value. To the extent that this plot is a 
straight line the distribution is normal. If the plot deviates greatly from a 
straight line the data does not have normal distribution. It may be seen that 
the plot in Figure 15 is a relatively straight line. A surer test is to use the 
log of probability and the log of the order value (Figure 16) . This is useful 
to test the values lying in the skirts of the distribution curve. Since there is 
always a small number of data points in the region of the lower adhesion values, 
the use of the log plot highlights any abnormality of these values. 
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Our analysis of these data have allowed us to conclude that the data is 
fairly normally distributed and can therefore be used to establish criteria for 
minimum threshold levels based on normal distribution. These criteria have been 
determined using an algorithm giving the probability that any number R of adhesion 
values will fall below the n lowest values. 

For the FAA we selected a test sample size of 20 tires and used the three 
lowest readings as the threshold criteria. Using the algorithm, we determined 
that, using values of 30, 33, and 36 for buff line adhesion and 20, 23, and 26 
for ply adhesion, the probability that a retreader having good tires would 
fail the test was about fifteen percent. The probability of failing a retest 
was about 2%. One the other hand the probability of detection of a sample of 
tires having a mean less than the threshold values increases very rapidly to 98% 
at a value of 1 standard deviation away from the threshold mean (Figure 17). 

Tread and ply adhesion values are a very good measure of tire production 
uniformity and can therefore be used as a monitor of quality during production in 
statistical QC devices such as control charts. The threshold values given rep- 
resent tires taken from a fleet in which a very small number of tire related 
incidents have occurred. They can therefore be considered as representing a safe 
population of tires. 
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VARIABLE: SIZE NUMBER: 2 DESCRIPTION; VARIABLE TYPE: FLOAT 

THERE HERE 706 OBSERVATIONS, HHICH INCLUDED 0 CASES OF HISSING DATA SELECTED FROM A TOTAL OF 706 OBSERVATIONS 


SUM OF OBSERVATIONS - 23559.43 

MEAN » 33.37030 

MAXIMUM - 56.16000 

STANDARD ERROR OF MEAN - 0.3228110 

COEFFICIENT OF SKEHNESS - 0.5019920 


SUM OF OBSERVATIONS SQUARED - 838052.0 

MEDIAN - 29.50750 

MINIMUM - 20.20000 

STANDARD DEVIATION - 8.577303 

COEFFICIENT OF VARIATION - 25.70341 


NUMBER OF OBSERVATIONS - 
MODE - 24.55000 

RANGE - 35.96000' 

VARIANCE » 73.57012 

KURTOSIS - 2.321392 


VALUE 

FREQUENCY 

PERCENTAGE 

CUMULATIVE 

PERCENTAGE 

20.20000 

1 

0.142 

0.142 

20.55000 

9 

1.275 

1.416 

22.55000 

50 

7.082 

8.499 

24.55000 

139 

19.688 

28.187 

26.66000 

17 

2.408 

30.595 

28.77000 

72 

10.198 

40.793 

28.90000 

65 

9.207 

50.000 

30.11500 

14 

1.983 

51.983 

36.11000 

94 

13.314 

65.297 

37.11500 

8 

1.133 

66.431 

40.14000 

129 

18.272 

84.703 

41.15000 

10 

1.416 

86.119 

44.13000 

18 

2.550 

88.669 

44.16000 

6 

1.133 

89.802 

46.14000 

2 

0.283 

90.085 

46.16000 

23 

3.258 

93.343 

47.18000 

1 

0.142 

93.484 

49.17000 

26 

3.683 

97.167 

50.20000 

3 

0.425 

97.592 

52.20000 

6 

0.850 

98.442 

56.16000 

11 

1.558 

100.000 
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Figure 2 
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VARIABLE: RL NUMBER: 3 DESCRIPTION: R LEVEL VARIABLE TYPE: FIICED 

THERE WERE 706 OBSERVATIONS, TffllCH INCLUDED 623 CASES OF MISSING DATA SELECTED FROM A TOTAL 706 OBSERVATIONS 

NUMBER OF OBSERVATIONS = 
MODE = 24.55000 

RANGE = 35.96000 

VARIANCE = 73.57012 

KURTOSIS = 2.321392 


VALUE 

FREQUENCY 

PERCENTAGE 

CUMULATIVE 

PERCENTAGE 

1 

19 

22.892 

22.892 

2 

13 

15.663 

38.554 

3 

17 

20.482 

59.036 

4 

16 

19.277 

78.313 

5 

11 

13.253 

91.566 

6 

3 

3.614 

95.181 

7 

1 

1.205 

96.386 

8 

1 

1.205 

97.590 

9 

1 

1.205 

98.795 

12 

1 

1.205 

100.000 
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SUM OF OBSERVATIONS =■ 269.00 

MEAN = 3.240964 

MAXIMUM - 12 

STANDARD ERROR OF MEAN = 0.2199233 
COEFFICIENT OF SKEWNESS » 1.395565 


SUM OF OBSERVATIONS SQUARED = 1201.000 

MEDIAN » 3.000000 

MINIMUM -= 1 

STAffl)ARD DEVIATION = 2.003597 

COEFFICIENT OF VARIATION = 61.82101 


Figure 3 
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VARIABLE: MFG NUMBER: A DESCRIPTION: MANUFACTURER VARIABLE TYPE: ALPHA 

THERE IffiRE 706 OBSERVATIONS, WHICH INCLUDED 430 CASES OF HISSING DATA SELECTED FROM A TOTAL OF 706 OBSERVATIONS 
MAXIMUM = TMS MINIMUM = AIR 


VALUE 

AIR 

BFG 

BRS 

DLP 

F 

GYR 

KC 

TMS 


FREQUENCY 

PERCENTAGE 

CUMULATIVE 

PERCENTAGE 

143 

51.812 

51.812 

34 

12.319 

64.130 

4 

1.449 

65.580 

7 

2.536 

68.116 

5 

1.812 

69.928 

36 

13.043 

82.971 

6 

2.174 

35.145 

41 

14.855 

100.000 
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Figure 4 
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VARIABLE: TMX NUMBER: 5 DESCRIPTION: MAX TREAD ADHSN VARIABLE TYPE: FIXED 

THERE WERE 706 OBSERVATIONS, IffllCH INCLUDED 563 CASES OF MISSING DATA SELECTED FROM A TOTAL OF 706 OBSERVATIONS 

NUMBER OF OBSERVATIONS = 
MODE = 80 

RANGE = 93 

VARIANCE = 408.1580 

KURTOSIS = 2.987374 


CU^aJLATIVE 


VALUE 


FREQUENCY 

PERCENTAGE 

PERCENTAGE 

50 - 

53 

5 

3.497 

3.497 

54 - 

57 

3 

2.098 

5.594 

58 - 

61 

6 

4.196 

9.790 

62 - 

65 

7 

4.895 

14.685 

66 - 

69 

12 

8.392 

23.077 

70 - 

73 

7 

4.895 

27.972 

74 - 

77 

11 

7.692 

35.664 

78 - 

81 

18 

12.587 

48.252 

82 - 

85 

10 

6.993 

55.245 

86 - 

89 

11 

7.692 

62.937 

90 - 

93 

7 

4.895 

67.832 

94 - 

97 

5 

3.497 

71.329 

98 - 

101 

10 

6.993 

78.322 

102 - 

105 

9 

6.294 

84.615 

106 - 

109 

5 

3.497 

88.112 

110 - 

113 

4 

2.797 

90.909 

114 - 

117 

2 

1.399 

92.308 

118 - 

121 

2 

1.399 

93.706 

122 - 

125 

1 

0.699 

94.406 

126 - 

129 

3 

2.098 

96.503 

130 - 

133 

1 

0.699 

97.203 

134 - 

137 

3 

2.098 

99.301 

138 - 

141 

0 

0.000 

99.301 

142 - 

145 

1 

0.699 

100.000 
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SUM OF OBSERVATIONS = 12289.00 

MEAN = 85.93706 

MAXIMUM = 143 

STANDARD ERROR OF MEAN = 1.609453 

COEFFICIENT OF SKEWIESS = 0.5453779 


SUM OF OBSERVATIONS SQUARED = 1114039. 

MEDIAN = 82.00000 

MINIMUM = 50 

STANDARD DEVIATION = 20.20292 

COEFFICIENT OF VARIATION = 23.50897 
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VARIABLE: TAV NUMBER: 7 DESCRIPTION: AVERAGE TREAD ADHSN VARIABLE TYPE: FIXED 

THERE HERE 706 OBSERVATIONS, WHICH INCLUDED 213 CASES OF MISSING DATA SELECTED FROM A TOTAL OF 


SUM OF OBSERVATIONS - 41810.00 

MEAN - 84.80730 

MAXIMUM » 149 

STANDARD ERROR OF MEAN - 0.9643130 

COEFFICIENT OF SKEWNESS - 0. 6466109E-01 


SUM OF OBSERVATIONS SQUARED - 3771346. 

MEDIAN = 85.00000 

MINIMUM “ 20 

STANDARD DEVIATION - 21.41122 

COEFFICIENT OF VARIATION - 25.24691 


706 OBSERVATIONS 

NUMBER OF OBSERVATIONS - 
MODE - 90 

RANGE - 129 

VARIANCE - 458.4405 

KURTOSIS » 2.747826 


CUMULATIVE 


VALUE 


FREQUENCY 

PERCENTAGE 

PERCENTAGE 

20 - 

25 

2 

0.406 

0.406 

26 - 

31 

0 

0.000 

0.406 

32 - 

37 

0 

0.000 

0.406 

38 - 

43 

7 

1.420 

1.826 

44 - 

49 

8 

1.623 

3.448 

50 - 

55 

26 

5.274 

8.722 

56 - 

61 

34 

6.897 

15.619 

62 - 

67 

34 

6.897 

22.515 

68 - 

73 

44 

8.925 

31.440 

74 - 

79 

38 

7.708 

39.148 

80 - 

85 

60 

12.170 

51.318 

86 - 

91 

52 

10.548 

61.866 

92 - 

97 

50 

10.142 

72.008 

98 - 

103 

38 

7.708 

79.716 

104 - 

109 

35 

7.099 

86.815 

110 - 

115 

29 

5.882 

92.698 

116 - 

121 

19 

3.854 

96.552 

122 - 

127 

3 

0.609 

97.160 

128 - 

133 

7 

1.420 

98.580 

134 - 

139 

4 

0.811 

99.391 

140 - 

145 

2 

0.406 

99.797 

146 - 

151 

1 

0.203 

100.000 
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Figure 6 
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VARIABLE; PMX NUMBER: 8 DESCRIPTION: MAX PLY ADHSN VARIABLE TYPE: FIXED 

THERE WERE 706 OBSERVATIONS, IffllCH INCLUDED 538 CASES OF MISSING DATA SELECTED FROM A TOTAL OF 


SUM OF OBSERVATIONS - 6691.000 

MEAN «= 39.82738 

MAXIMUM » 73 

STANDARD ERROR OF MEAN - 0.7058958 

COEFFICIENT OF SKEWNESS - 1.399312 


SUM OF OBSERVATIONS SQUARED - 280465,0 

MEDIAN - 38.00000 

MINIMUM - 26 

STANDARD DEVIATION - 9.149456 

COEFFICIENT OF VARIATION ■= 22.97270 


706 OBSERVATIONS 

NUMBER OF OBSERVATIONS 
MODE = 39 

RANGE = 47 

VARIANCE = 83.71254 
KDRTOSIS - 5.195022 


VALUE 


FREQUENCY 

PERCENTAGE 

CUMULATIVE 

PERCENTAGE 

26 - 

27 

5 

2.976 

2.976 

28 - 

29 

6 

3.571 

6.548 

30 - 

31 

12 

7.143 

13.690 

32 - 

33 

15 

8.929 

22.619 

34 - 

35 

20 

11.905 

34.524 

36 - 

37 

20 

11. 905 

46.429 

38 - 

39 

23 

13.690 

60.119 

40 - 

41 

13 

7.738 

67.857 

42 - 

43 

14 

8.333 

76.190 

44 - 

45 

9 

5.357 

81.548 

46 - 

47 

8 

4.762 

86.310 

48 - 

49 

3 

1.786 

88.095 

50 - 

51 

2 

1.190 

89.286 

52 - 

53 

4 

2.381 

91.667 

54 - 

55 

0 

0.000 

91.667 

56 - 

57 

3 

1.786 

93.452 

58 - 

59 

3 

1.786 

95.238 

60 - 

61 

1 

0.595 

95.833 

62 - 

63 

1 

0.595 

96.429 

64 - 

65 

1 

0.595 

97.024 

66 - 

67 

1 

0.595 

97.619 

68 - 

69 

1 

0.595 

98.214 

70 - 

71 

2 

1.190 

99.405 

72 - 

73 

1 

0.595 

100.000 
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Figure 7 
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VARIABLE: PMN NUMBER: 9 DESCRIPTION; MIN PLY ADHSN VARIABLE TYPE: FIXED 

THERE WERE 706 OBSERVATIONS, WHICH INCLUDED 681 CASES OF HISSING DATA SELECTED FROM A TOTAL OF 


SUM OF OBSERVATIONS - 711.0000 

MEAN - 28.4A000 

MAXIMUM - 44 

STANDARD ERROR OF MEAN - 1.605283 

COEFFICIENT OF SKEWNESS - 0.3450782 


SUM OF OBSERVATIONS SQUARED » 21767.00 

MEDIAN - 27.00000 

MINIMUM - 14 

STANDARD DEVIATION - 8.026415 

COEFFICIENT OF VARIATION - 28.22227 


706 OBSERVATIONS 

NUMBER OF OBSERVATIONS = 
MODE » 20 

RANGE - 30 

VARIANCE - 64.42333 

KURTOSIS = 2.162221 


VALUE 

FREQUENCY 

PERCOTTAGE 

CUMULATIVE 

PERCENTAGE 

14 

1 

4.000 

4.000 

19 

1 

4.000 

8.000 

20 

3 

12.000 

20.000 

22 

2 

8.000 

28.000 

24 

2 

8.000 

36.000 

25 

1 

4.000 

40.000 

26 

2 

8.000 

48.000 

27 

2 

8.000 

56.000 

29 

1 

4.000 

60.000 

30 

2 

8.000 

68.000 

31 

1 

4.000 

72.000 

36 

1 

4.000 

76.000 

37 

1 

4.000 

80.000 

38 

2 

8.000 

88.000 

40 

1 

4.000 

92.000 

42 

1 

4.000 

96.000 

44 

1 

4.000 

100.000 
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Figure 8 
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VARIABLE: PAV NUMBER: 10 DESCRIPTION: AVERAGE PLY ADUSN VARIABLE TYPE: 

THERE WERE 706 OBSERVATIONS, WHICH INCLUDED 188 CASES OF MISSING DATA SELECTED FROM 


SUM OF OBSERVATIONS - 28109.00 

MEAN “ 54.26448 

MAXIMUM = 114 

STANDARD ERROR OF MEAN - 0.5346928 

COEFFICIENT OF SKEWNESS - 0.6843333 


SUM OF OBSERVATIONS SQUARED » 1601885. 

MEDIAN - 54.00000 

MINIMUM - 17 

STANDARD DEVIATION » 12.16940 

COEFFICIENT OF VARIATION » 22.42609 


VALUE 
17 - 

20 

FREQUENCY 

1 

PERCENTAGE 

0.193 

CUMULATIVE 

PERCENTAGE 

0.193 

21 - 

24 

0 

0.000 

0.193 

25 - 

28 

3 

0.579 

0.772 

29 - 

32 

9 

1.737 

2.510 

33 - 

36 

19 

3.668 

6.178 

37 - 

40 

37 

7.143 

13.320 

41 - 

44 

36 

6.950 

20.270 

45 - 

48 

41 

7.915 

28.185 

49 - 

52 

77 

14.865 

43.050 

53 - 

56 

92 

17.761 

60.811 

57 - 

60 

62 

11.969 

72.780 

61 - 

64 

52 

10.039 

82.819 

65 - 

68 

42 

8.108 

90.927 

69 - 

72 

23 

4.440 

95.367 

73 - 

76 

5 

0.965 

96.332 

77 - 

80 

5 

0.965 

97.297 

81 - 

84 

4 

0.772 

98.069 

85 - 

88 

3 

0.579 

98.649 

89 - 

92 

1 

0.193 

98.842 

93 - 

96 

2 

0.386 

99.228 

97 - 

100 

1 

0.193 

99.421 

101 - 

104 

0 

0.000 

99.421 

105 - 

108 

2 

0.386 

99.807 

109 - 

112 

0 

0.000 

99.807 

113 - 

114 

1 

0.193 

100.000 
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Figure 9 


FIXED 

A TOTAL OF 

706 OBSERVATIONS 

NUMBER OF OBSERVATIONS “ 
MODE - 53 

RANGE “ 97 

VARIANCE - 148.0943 

KURTOSIS « 5.372330 
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VARIABLE: TENS NUMBER: 10 DESCRIPTION: TENSILE VARIABLE TYPE: FIXED 

THERE WERE 706 OBSERVATIONS, WHICH INCLUDED 134 CASES OF MISSING DATA SELECTED FROM A TOTAL OF 


SUM OF OBSERVATIONS “ 1538885. 

MEAN - 2690.358 

MAXIMUM - 3697 

STANDARD ERROR OF MEAN - 13.54925 

COEFFICIENT OF SKEWNESS - 0.3102679 


SUM OF OBSERVATIONS SQUARED = 
MED UN - 2661.000 

MINIMUM - 1495 

STANDARD DEVUTION = 324.0510 
COEFFICIENT OF VARUTION - 12 


0.4200112E+10 OBSERVATIONS 

NUMBER OF OBSERVATIONS » 
MODE - 2660 

RANGE » 2202 

04490 VARIANCE •= 105009.1 

KURTOSIS - 3.827942 


VALUE 


FREQUENCY 

PERCENTAGE 

CUMULATIVE 

PERCENTAGE 

1495 

- 

1583 

1 

0.175 

0.175 

1584 

- 

1672 

1 

0.175 

0.350 

1673 

- 

1761 

1 

0.175 

0.524 

1762 

- 

1850 

1 

0.175 

0.699 

1851 

- 

1939 

2 

0.350 

1.049 

1940 

- 

2028 

5 

0.874 

1.923 

2029 

- 

2117 

5 

0.874 

2.797 

2118 

- 

2206 

14 

2.448 

5.245 

2207 

- 

2295 

17 

2.972 

8.217 

2296 

- 

2384 

31 

5.420 

13.636 

2385 

- 

2473 

58 

10.140 

23.776 

2474 

- 

2562 

65 

11.364 

35.140 

2563 

- 

2651 

75 

13.112 

48.252 

2652 

- 

2740 

85 

14.860 

63.112 

2741 

- 

2829 

50 

8.741 

71.853 

2830 

- 

2918 

40 

6.993 

78.846 

2919 

- 

3007 

24 

4.196 

83.042 

3008 

- 

3096 

29 

5.070 

88.112 

3097 

- 

3185 

30 

5.245 

93.357 

3186 

- 

3274 

11 

1.923 

95.280 

3275 

- 

3363 

10 

1.748 

97.028 

3364 

- 

3452 

3 

0.524 

97.552 

3453 

- 

3541 

5 

0.874 

98.427 

3542 

- 

3630 

4 

0.699 

99.126 

3631 

- 

3697 

5 

0.874 

100.000 
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Figure 10 
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VARIABLE: ELONG NUMBER: 


12 DESCRIPTION: ELONGATION 


VARIABLE TYPE: FIXED 


THERE WERE 706 OBSERVATIONS, WHICH INCLUDED 133 CASES OP MISSING DATA SELECTED FROM A TOTAL OF 


SUM OF OBSERVATIONS » 281707.0 

MEAN - 491.6353 

MAXIMUM =■ 670 

STANDARD ERROR OF MEAN 
COEFFICIENT OF SKEWNESS 


2.767603 

-0.1577621 


SUM OF OBSERVATIONS SQUARED 
MEDIAN - 497.0000 , 

MINIMUM =■ 274 

STANDARD DEVIATION - 66.24928 
COEFFICIENT OF VARIATION = 13.47529 


0.1410076E+09 


706 OBSERVATIONS 

NUMBER OF OBSERVATIONS 
MODE = 530 

RANGE = 396 

VARIANCE = 4388.966 

KURTOSIS = 2.601573 


CUMULATIVE 


VALUE 


FREQUENCY 

PERCENTAGE 

PERCENTAGE 

274 - 

289 

1 

0.175 

0.175 

290 - 

305 

1 

0.175 

0.349 

306 - 

321 

0 

0.000 

0.349 

322 - 

337 

1 

0.175 

0.524 

338 - 

353 

4 

0.698 

1.222 

354 - 

369 

15 

2.618 

3.839 

370 - 

385 

10 

1.745 

5.585 

386 - 

401 

25 

4.363 

9.948 

402 - 

417 

33 

5.759 

15.707 

418 - 

433 

44 

7.679 

23.386 

434 - 

449 

22 

3.839 

27.225 

450 - 

465 

35 

6.108 

33.333 

466 - 

481 

49 

8.551 

41.885 

482 - 

497 

48 

8.377 

50.262 

498 - 

513 

59 

10.297 

60.558 

514 - 

529 

47 

8.202 

68.761 

530 - 

545 

52 

9.075 

77.836 

546 - 

561 

46 

8.028 

85.864 

562 - 

577 

23 

4.014 

89.878 

578 - 

593 

31 

5.410 

95.288 

594 - 

609 

9 

1.571 

96.859 

610 - 

625 

9 

1.571 

98.429 

626 - 

641 

6 

1.047 

99.476 

642 - 

657 

0 

0.000 

99.476 

658 - 

670 

3 

0.524 

100.000 


573 
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VARIABLE: DUR NUMBER! 17 DESCRIPTION: DUROMETER VARIABLE TYPE: FLOAT 

THERE WERE 706 OBSERVATIONS. WHICH INCLUDED 276 CASES OF MISSING DATA SELECTED FROM A TOTAL OF 


SUM OF OBSERVATIONS - 26392.00 

MEAN - 61. 37 67 A 

MAXIMUM - 69.00000 

STANDARD ERROR OF MEAN - 0.1202394 

COEFFICIENT OF SKEWNESS - -0.4979397 


SUM OF OBSERVATIONS SQUARED - 1622522. 

MEDIAN - 61.00000 

MINIMUM - 50.00000 

STANDARD DEVIATION - 2.493337 

COEFFICIENT OF VARIATION - 4.062348 


706 OBSERVATIONS 

NUMBER OF OBSERVATIONS - 
MODE » 61.00000 

RANGE - 19.00000 
VARIANCE - 6.216729 

KURTOSIS - 5.129465 


VALUE 

FREQUENCY 

PERCENTAGE 

CUMULATIVE 

PERCENTAGE 

50.00000 

2 

0.465 

0.465 

54.00000 

1 

0.233 

0.698 

55.00000 

9 

2.093 

2.791 

56.00000 

2 

0.465 

3.256 

57.00000 

13 

3.023 

6.279 

58.00000 

14 

3.256 

9.535 

59.00000 

28 

6.512 

16.047 

60.00000 

70 

16.279 

32.326 

61.00000 

96 

22.326 

54.651 

62.00000 

56 

13.023 

67.674 

63.00000 

63 

14.651 

82.326 

64.00000 

38 

8.837 

91.163 

65.00000 

24 

5.581 

96.744 

66.00000 

8 

1.860 

98.605 

67.00000 

2 

0.465 

99.070 

68.00000 

1 

0.233 

99.302 

69.00000 

3 

0.698 

100.000 


430 
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Figure 12 
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VARIABLE: 


SG 


innCBER: 


18 


DESCRIPTION: TREAD SPECIFIC GRAV l.()() VARIABLE TYPE: 


float 


THERE HERE 706 OBSERVATIONS, WHICH INCLUDED 279 CASES OF MISSING DATA SELECTED FROM A TOTAL OF 


SUM OF OBSERVATIONS - 5376.000 

MEAN - 12.59016 

MAXIMUM •= 16.00000 

STANDARD ERROR OF MEAN ^ 
COEFFICIENT OF SKEWNESS 


0.5154423E-01 

-0.2478130 


68168.00 


SUM OF OBSERVATIONS SQUARED 
MEDIAN = 13.00000 

MINIMUM » 8.000000 

STANDARD DEVIATION = 1.065109 

COEFFICIENT OF VARIATION » 8.459850 


706 OBSERVATIONS 

NUMBER OF OBSERVATIONS 
MODE = 12.00000 

RANGE = 8.000000 

VARIANCE = 1.134457 

KURTOSIS • 4.195873 


VALUE 

FREQUENCY 

PERCENTAGE 

CUMULATIVE 

PERCENTAGE 

8.000000 

1 

0.234 

0.234 

9.000000 

3 

0.703 

0.937 

10.00000 

7 

1.639 

2.576 

11.00000 

34 

7.963 

10.539 

12.00000 

161 

37.705 

48.244 

13.00000 

143 

33.489 

01.733 

14.00000 

65 

15.222 

96.956 

15.00000 

12 

2.810 

99.766 

16.00000 

1 

0.234 

100.000 


427 




10.0 20.0 30.0 40.0 

I 

IX 

IXX 

IXXXXXXXX 

IXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

IXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

IXXXXXXXXXXXiCXXX 

IXXX 

I 

H 1 h H h— 

10.0 20.0 30.0 40.0 


50.0 


50.0 


Figure 13 
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50% 



-3 -2 -1 MEAN OF +1 +2 +3 +4 

ALL VALUES 


Figure 14 


PLOT 

OF VARIABLE; 1 (llORIZ.) VS VARIABLE: 3 

(VERT.) 

1.000 

+ 

3B 



I 

CA 



I 

4CD 



I 

AC 9 


0.8000 

+ 

3C3 



I 

2CA 



I 

ICB 



I 

2CC 


0.6000 

+ 

ACB 



I 

93 



I 

2CD3 



I 

CB 


0.4000 

+ 

6D 



I 

3D6 



I 

79 



I 

CC6 


0.2000 

+ 

3DC1 



I 

79 



I 

6DC6 



I BD6 

0.1863E-07 +71 

I 
I 
I 

- 0.2000 + 


201.5 451.5 701.5 

326.5 576.5 826.5 


Figure 15 
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PROBABILITY 
THAT A BAD 
LOT WILL BE 
DETECTED 


PLOT OF VARIABLES VLOGl (HORIZ.) VS VARIABLE: VLOG3 (VERT.) 


I 

0,0000 +1 
I 
I 
I 

-O.AOOO + 

I 

I 

I 

-0.8000 + 

I 

I 

I 

- 1.200 + 

I 

I 

I 

-1.600 + 

I 

I 

I 

- 2.000 + 

I 

I 

I 

-2.400 + 

I 

-0.8000E-01 


1 


H — 4 : 1 h 

YF 

** 

3*1 

B* 

18 

8E 

BJ 

IB 

A 

66 

9B6 

5 

5562 

13 

12 

131 


112 


1 1 


0,7200 


1.520 


2.320 


3.120 


Figure 16 
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HEAT GENERATION IN AIRCRAFT TIRES 

Samuel K. Clark 
University of Michigan 

ABSTRACT 

A method has been developed for calculating the internal temperature distribu- 
tion in an aircraft tire while free rolling under load. The method uses an approxi- 
mate stress analysis of each point in the tire as it rolls through the contact 
patch. From this stress change, the mechanical work done on each volume element 
may be obtained and converted into a heat release rate through a knowledge of 
material characteristics- The tire cross-section is then considered as a body with 
internal heat generation, and the diffusion equation is solved numerically with 
appropriate boundary conditions of the wheel and runway surface. Comparison with 
data obtained with burled thermocouples in tires shows good agreement. The 
data presented in this paper were excerpted from reference 1. 
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The origin of aircraft tire heating lies in the hysteretic loss characteristics 
of pol 3 nnerlc materials used for both the nylon textile cord reinforcement and for 
the rubber components in an aircraft tire. Figure 1 shows a typical stress- 
strain curve between the two end points of the cyclic stress states associated with 
rolling a material element through the contact patch of the tire and out again. 



€ 


Figure 1 
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The stress excursion is presumed to occur between two points A and B located 
180“ apart on the rolling tire. One extreme is the upper point A where only in- 
flation pressures act to produce an axisymmetric stress state of the tire. The 
other extreme point is the center of the contact patch, where again an axisymmetric 
solution is used to obtain the stress state in the tire deformed against a flat 
surface (fig. 2). 
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Figure 3 Illustrates the assumed conditions in the vicinity of the contact 
patch of the tire. A total deflection A is presumed to take place between the 
outer surface of the tire crown and the ground. Inside this contact patch region 
the tire is presumed to be in complete contact with the flat runway surface. 



(a) Assumed geometry at point A 
of figure 2. 


(b) Assumed geometry at point B 
of figure 2. 


Figure 3 
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Figure 4 shows the assumed geometry of the tire at the two extreme condi- 
tions A and B illustrated In figures 2 and 3. In the case of axlsymmetric infla- 
tion loads the geometric midline of the tire is represented by a torus as shown. 
For that portion of the tire in the center of the contact patch, the shape is 
considered to be flat in the contact region and made up of a single circular arc 
outside the contact region. The stress solutions are both obtained from axisym- 
metrlc considerations. 




197 


The tire is divided into five segments for purposes of both stress and thermal 
analysis (fig. 5). Segment I is that portion of the tire carcass lying in the contact 
patch. It has plane orthotropic elastic properties. Segment II is isotropic 
tread material in contact with the runway surface, and hence is acted upon by 
normal pressures as well as by membrane and bending strains. Segment III is that 
part of the side wall acted on by internal pressure, and also has orthotropic elas- 
tic properties. Segment IV is the sidewall cover of the tire and is considered to 
be isotropic in its material properties. Segment V is the bead region which is 
subject to a special description. 


Edge of 



Figure 5 
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Figure 6 shows the method of splitting the tire cross section into various 
rectangular segments for purposes of thermal analysis. In general the decomposi- 
tion process follows the natural boundary lines of the tire between the orthotropie 
and isotropic regions. 
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Figure 7 is the general heat diffusion equation which is being modeled on a 
discrete basis in the numerical calculation described here. However, the elements 
which are in contact with the outer surface or with the bead area of the tire have 
additional terms associated with either thermal conduction into the metallic wheel 
or convective heat transfer loss to the air. 


K e 


p 


8t 


q 


where 

K = thermal conductivity 
p = density 

o 

q = heat generation rate 
0 = temperature 


Figure 7 
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The eqtiation for convective heat transfer, temperature buildup and internal 
heat generation is discretized into the form shovm in figure 8 for numerical 
calculation throughout the tire cross section. Note that due to the higher ro- 
tational speeds of the tire compared to thermal diffusion times, the temperature 
distribution is completely axisymmetrlc. 


^ “k,K+1,£ 

“k, £-!,£ ^ ®K, 

^ ®K, £+1“®K, 

with 

KA 

l,m,n d 

K = thermal conductivity 

A = area of contact between subscripted elements 

d = distance between centroids of subscripted 
elements 

NOTATION FOR THERMAL ELEMENTS 
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The strain energy density must be related to the hysteretlc loss, which results 
in internal heat generation. Due to the lack of appropriate material constants, 
this is now being done on a provisional basis using the loss tangent tan 6 , a quan- 
tity usually associated with uniaxial sinusoidal loading. Due to the linear rela- 
tionship between rate of heat generation and frequency, the rate of heat generation 
is linearly proportional to aircraft velocity (fig. 9). 


~ n tan5 , . 3 , 

q = 0.01 cal/cmVsec 

whare 

Tq = outside radius of tire, in. 
A = tire deflection, in. 

Vq = aircraft velocity, ft/sec 
U = elastic energy, in.-lb/in.3 

Figure 9 
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In data shown In subsequent figures, calculations are carried out comparing 
predictions of this program with measured temperature rise data. These predictions 
are based on best estimates of material properties taken from the literature. 

(ref. 2). These are given in figure 10. 


MATERIAL CHARACTERISTICS OF 
22 X 5.5 and 40 X 14 AIRCRAFT TIRES 


Region^ 

Young ' s 
modulus 
of rubber , 
E, psi 

Shear 
modulus 
of rubber , 
G, psi 

*’tan6o, 

psi 

Cord 
angle, 
a, deg 

Thermal 

conductivity, 

K, 

cal-cm/“C-cm -sec 

Density, 
P. deg 

Specific 

heat, 

Cp, 

Tread (II) 


335 

0.15 


5 X lo"^ 

1.0 

0.5 

Carcass 
(I, III) 

(c) 


0.15 

(d) 

5 X 10“^ 

1.0 

0.5 

Sidewall 
rubber (IV) 


335 

0.15 



1.0 

0.5 

Bead (V) 

15 X 10^ 


0.03 






^See fig. 5. 

^Values of tan6g are given for 25°C. For higher temperatures, tan6 is calculated 
by the approximate expression 

tan6 = tan6^(e"*°^^®^ 

where A0 is the temperature rise in °C above ambient. 

*^See reference 1. For shear modulus G of rubber use 335 psi. 

^Calculated by cosine law using a = 35“ crown angle. 


Figure 10 
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Figure 12 shows comparison between temperature rise and time as measured by 
thermocouples, along with calculations for similar locations on the tire using the 
computer program previously described. Stations A through F correspond to those 
shown in figure 11. 


TIRE:22x5.5-8PR 
SPEED:20MPH 
FREE ROULING 
UM TEST(REF. 1) 


O INSIDE 
A MIOLINE 
□ OUTSIDE 
— CALCULATED 


SURFACE; 120 IN. DIA.DRUM 
Fz ; 4350/4350 = 100% RATED 
8z : 1.341 / 1. 344 = 99.8 % rated 
P o: 125/115 =109% RATED 



Figure 12 
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TEMPERATURE RISE,°F 


Figure 13 shows similar temperature vs. time information for a different form 
of this 22 X 5.5 tire. In this case the tire is a 12 PR tire and requires a dif- 
ferent geometric description due to the thicker carcass. This in turn results in 
considerably higher temperatures being generated as this tire is run under rated 
load conditions for approximately 150 seconds. Again, measured data is compared 
with calculations from the computer program. 


TIRE i 22x5.5 12 PR O INSIDE SURFACE : 120 IN STEEL DRUM 

SPEED:20MPH ^ MIDLINE Fz 7J00/7I00 100% RATED 

FREE ROLLING □ OUTSIDE 8z = 1-304/1.312 99% RATED 

UMTEST(REF. 1) CALCULATED Po = 209/190 110%RATED 



TIME, SEC 


Figure 13 
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TEMPERATURE RISE 


Figure 14 shows a cross section of a 40 x 14 28 PR tire used as a model of a 
large aircraft tire. This is a common commercial tire, and in this case a larger 
number of discrete sections is chosen for temperature calculation purposes. Again, 
thermocouple locations are shown by the black dots. Extensive experimental work 
has been done on this tire in order to compare it with the calculation. 


O 


TIRE;40xl4-22PR 
SPEED; 20 MPH 
FREE ROLLING 
UM TEST(REF. 1) 


O INSIDE 
A MIDLINE 
a OUTSIDE 
— CALCULATED 


SURFACE: 120 IN. DIA. DRUM 
Fz ; 23,500 / 25,000 = 94 % RATED 
8z : 3.378 / 3.378 = 100 % RATED 
Po ; 160 / 155 « 103 % RATED 



TIME, SEC 


Figure 14 
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Figure 15 shows a. comparison between a 22 ply rated version of the tire shown 
in figure 14 and calculations carried out using the program described previously. 
These calculations generally agree will with measured data. 



Dimensions in inches 

Figure 15 
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Figure 16 shows temperature rise data up to approximately 150 seconds of run- 
ning time under essentially rated conditions for the 28 ply rated version of this 
tire. In this case calculations agree well with experimental measurements, although 
here the thicker section tire also requires a separate geometric description. 


TIRE; 40XI4-28PR O INSIDE SURFACE; 120 IN.DIA. DRUM 

SPEED: 20 MPH A MIDLINE Fz . 30,800/ 33,500 = 92 % RATED 

FREE ROLLING □ OUTSIDE : 3.459/3.484 = 99 % RATED 

UM TEST (REF. 1) — CALCULATED Pq ; 207/200= 104 % RATED 



50 100 150 0 50 100 150 c 50 100 150 c 50 100 150 0 50 100 150 0 50 100 150 


TIME, SEC 
Figure 16 
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RESULTS FROM RECENT NASA TIRE THERMAL STUDIES 


John L. McCarty 
NASA Langley Research Center 


ABSTRACT 

This paper describes the testing technique and some results from an experimental 
study to determine tire temperature profiles to aid in defining the strength and 
fatigue limitations of the tire carcass structure. This effort is part of a program 
to explore analytically and through experiment the temperature distribution in an 
aircraft tire during free roll and braked and yawed rolling conditions. The analytical 
effort, together with a comparison with the experimental results, is discussed in 
the paper by Clark in this publication (ref. 1). 


LOCATION OF THERMOCOUPLES IN TIRE THERMAL STUDY 


Figure 1 shows the approximate locations of the eighteen thermocouples installed 
within the carcass of several tires to support the thermal study. As noted in the 
figure, the thermocouples were mounted on the inner and outer walls of the tires and 
along an approximate midline at six radial stations. Most of the thermocouple 
installation was accomplished by implanting the sensors in holes drilled into the 
tire carcasses after the carcasses had been buffed prior to retreading. 



Figure 1 


OilSimi” f§ 

OF QUALITY 


PHOTOGRAPH OF AN INSTRUMENTED TEST TIRE 

Figure 2 Is a photograph of a tire ready for testing, and shows how the outputs 
from the thermocouples were routed through a slip ring en route to a data logger 
located in the instrument section of the test vehicle pictured in figure 3. The 
tires for these tests were size 22 x 5.5, 12-ply rating, type VII aircraft tires. 



Figure 2 
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INSTRUMENTED TIRE TEST VEHICLE 


ORIGINAL PAGE IS 
OF POOR QUALITY 


Figure 3 is a photograph of the powered ground test vehicle used in this jnves- 
tigatior This vehicle consists of a truck which has been to accommodate 

a tire test fixture and supporting equipment. For the tests described here, the 
tires were vertically loaded to 4000 pounds and driven at 20 mph, except for one 
series where speeds were increased to 50 mph. over a 
measurements recorded every 10 seconds. Data were taken 7’'^. 

yawed rolling, and at fixed slip ratios to simulate braked-rolling conditions. 



Figure 3 
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TYPICAL TEMPERATURES THROUGH THE TIRE CARCASS 
DURING FREE ROLLING 


Figure 4 presents the time history of temperature buildup In six of the 
thermocouples as the tire In the free-rolling mode made two and a half passes of 
about 6500 feet each down an asphalt runway. As noted on the figure, the test tire 
was raised from the pavement surface following each pass while the test vehicle 
was being turned around for the return trip. The thermocouples were selected to 
Illustrate the temperature rise through the tire carcass In the shoulder and bead 
areas. The figure clearly demonstrates that In both areas the rate of temperature 
buildup, and hence the magnitude of the temperature. Is greatest along the Interior 
wall, decreases through the carcass, and Is lowest along the exterior wall. As 
would be expected because of tire flexure, the temperatures In the shoulder area 
are greater than the corresponding temperatures In the bead area. 


TIRE; 22x5.5, 12 PLY TIRE DEFL. = 25% V^= 20 MPH 



Figure 4 
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TYPICAL INTERIOR TIRE CARCASS TEIIPERATURES DURING 

FREE ROLL 

Temperature data from the six thermocouples installed along the tire internal 
wall are presented in figure 5 as a function of the free-roll distance. It is 
interesting to note that the maximum temperature was measured not in the center of 
the tread nor in the shoulder area, but between the two. The temperature in the 
bead area is again shown to be well below the other measured tire internal tempera 
tures . 


22 X 5. 5, 12-PLY 


TIRE DEFL = 


V = 20 MPH 

g 


TIRE OFF GROUND- 


DISTANCE, FT 


Figure 5 




EFFECT OF TIRE DEFLECTION 


Figure 6 Illustrates the effect of tire deflection 6 on the distribution 
of temperature within the tire carcass after the tire had free— rolled for a distance 
of 5000 feet. The different deflections in this case were obtained by changing the 
tire inflation pressure while maintaining the 4000-pound vertical loading. Although 
data are lacking in the tread region, similar tests on other tires showed little 
effect of tire deflection on temperatures in that area. However, the figure shows 
a pronounced effect on temperatures in the sidewall region. The higher temperatures 
associated with the greater deflections are attributed to the more severe sidewall 
flexing which Increases with tire deflection- Similar trends were observed when 
the tire deflection was changed by varying the vertical load while maintaining a 
fixed Inflation pressure. 


FREE ROLL d 



5000' . 'p:;! 



C (°F) 

! J <27 

( <80) 

27 to 38- 

( 80 to 100) 

■i 38. to 49 

(100 to 120): 


(120 to 140): 

■1 60 to 71 

(140 to 160) 

IH 71 to 82 

(160 to 180): 


( >180) 




Figure 6 
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EFFECT OF YAW ANGLE 


To examine the effect of yaw angle ijj on tire temperature distribution, a series 
of free-rolling tests was conducted where the tire was yawed to 3° and 6°. Because 
only one side of the tire was instrumented, two passes were required at each yaw 
angle with the tire yawed first in one direction (+3°, for example) and then in 
the other direction (-3°, for example) to complete the temperature picture. The 
results are illustrated in figure 7 which presents the tire temperature profiles 
after having traveled 7000 feet. The figure shows a considerable temperature build- 
up with yaw angle in the tread area and in the "downwind" sidewall of the tire with 
an accompanying cooling effect on the "upwind" side. The reduced temperature on the 
"upwind" side is apparently due to the stress relieving which occurs in the "upwind" 
sidewall when the shoulder area on that side is forced to roll partially into the 
footprint. Conversely, the greater tire distortion on the "downwind" side as the 
yaw angle is Increased would account for those higher temperatures. The higher 
temperatures in the tread area can be explained by the increased scrubbing action 
in the tire footprint that is associated with increasing yaw angle. 


6 = 30%; d 






mj 


’■"'1 125) 


52,; to -6,6 

<125 to 150) 

■i. 

66 to 79 

<150 to 175) 

■■ 

79 to ‘93 

(175 to 200) 

■P; 

93 to 107 

(200 to 225) 

■■ 

107 to 121 

(225 to 250) 

■1 

>121 

( >250) 


Figure 7 
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EFFECT OF BRAKING 


Driving the test tire at fixed slip ratios provides a means for evaluating the 
effect of braking on the temperature buildup in the tire carcass. Figure 8 presents 
the temperature profiles after the tire had traveled 7000 feet first in free roll 
and subsequently at slip ratios of 5 and 10 percent. The figure shows that the 
influence of braking is essentially limited to the tread area. Minimal changes 
are noted in the temperature profiles of the sidewall regions. 




°C ('^F) 


C~i <49 

( <120) 

49 to 60 

(120 to 140) 

H 60 to 71 

(140 to 160) 

■■ 71 to 82 

(160 to 180) 

■i 82 to 93 

(180 to 200) 

■1 93 to 104 

(200 to 220) 

HI > 3^04 

{ > 220) 


Figure 8 
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EFFECT OF SPEED 


Over the speed range of these tests, the influence of wheel speed on the tire 
temperature distribution is almost insignificant. As observed in figure 9, the 
sidewall temperature after free-rolling 7000 feet is essentially unaffected and only 
moderate Increases can be seen in the tread temperature as the speed is increased. 
The rise in tread temperature at higher speeds is perhaps the result of the higher 
frequency squirming action which has been observed in the tire footprint. 


FREE ROLL 6 = 30%; d = 7000' 





□ <49 

( < 120) 

■n 

82 to 93 

(180 to 200) 

i 1 49 to 60 

(120 to 140) 

■i 

93 to 104 

(200 to 220) 

60 to 71 

<140 to 160) 

■i 

>104 

( > 220) 

71 to 82 

(160 to 180) 





Figure 9 
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TIRE EQUILIBRIUM TEMPERATURES 


The question arose during the course of this investigation as to whether the 
tire carcass would assume an equilibrium temperature profile if permitted to free 
roll for an indefinite period. To answer the question, many passes were made in 
both directions along a 7000-foot runway at 20,mph with a very brief turn-around 
time between passes to minimize any tire cooling. It became apparent after the 
tire had traveled approximately 65,000 feet that each thermocouple was approaching 
an equilibrium value. The final temperature profile, obtained after the tire had 
traveled in excess of 80,000 feet, is presented in figure 10 and shows temperatures 
approaching 300° F near the shoulder area. Note again that, particularly in the 
sidewall region, the temperatures along the inner wall are considerably greater than 
those near the outer surface. Also note that the centerline tread temperatures are 
somewhat lower than other temperatures in the tread. 
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CONCLUDING COMMENTS 


The results from a program to study tire temperature profiles to aid in 
defining the strength and fatigue limitations of an aircraft tire carcass struc- 
ture suggest the following concluding comments. 


• A testing technique has been developed for successfully measuring the 

thermal characteristics of a rolling tire. 

• Tire shoulder and sidewall temperatures increase with increasing tire 

deflection due to either overload or underinflation. 

• Tire tread and "downwind" sidewall temperatures Increase with yaw 

angle. 

• Large temperature increases are noted in the tread region during 

increased braking. 

• Changes in ground speed (over the test range to 50 mph) produce 

insignificant effects on the temperature buildup in the tire carcass. 
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SUMMARY OF SESSION DISCUSSIONS 


John A. Tanner 
INTRODUCTION 

General discussion periods were held at the end of each session to allow the 
participants to describe their own findings or to define problem areas In which they 
need help. In addition, a panel discussion was held among the session chairmen at 
the end of the workshop to review the status of tire modeling. These discussions 
Identified two levels of tire modeling that were considered significant. One level 
was that required by the vehicle dynamicist who is concerned with the vehicle and 
tire combination and the response of the vehicle to braking and steering Inputs, 
shimmy, and pavement roughness. The other level of modeling was that required by 
the tire designer who is concerned with detailed information on the stresses and 
strains associated with tire loading, carcass temperature profiles, and properties 
of the tire constituents. The following paragraphs will assess the status of tire 
modeling from the points of view of the vehicle dynamicist and the tire designer. 

STATUS OF TIRE MODELING - VEHICLE DYNAMICIST POINT OF VIEW 

The vehicle dynamicist is typically interested in the interaction between tire 
behavior and the performance of vehicular suspension systems or aircraft landing 
gear systems. For any simulation study, the dynamicist is generally seeking such 
tire response characteristics as surface interface friction, spring rates, 
hysteresis loss properties, and relaxation lengths. Ideally, these properties 
can only be obtained from experimental measurements, and this process may require 
considerable computer storage if the dynamicist is interested in simulations in- 
volving many different tire sizes or designs. Hoxjever the use of semiemplrlcal 
methods for predicting tire properties has met with some success. References 1 
and 2 present examples for estimating mechanical properties of various aircraft 
tires . 

The finite element method is generally considered to be too costly for this 
type of tire modeling, and the vehicle dynamicist typically resorts to much simpler 
models of the tire. The millions of well-behaved suspension systems and aircraft 
landing gear systems which provide good braking and cornering response and effec- 
tively eliminate shimmy and other potentially hazardous dynamic vibrations demon- 
strate the effectiveness of these simplified tire models. 

From the point of view of the vehicle d 3 mamicist, there is room for in- 
creased cooperation between the vehicle dynamicists in the automobile industry 
and the landing gear dynamicists in the aircraft industry. Although both groups 
of dynamicists work on similar tire problems, they tend to use different nomen- 
clature. This fact suggests that the two groups generally do not communicate 
with each other or read each other's publications, and such a failure in communi- 
cation can lead to a duplication of effort. 

STATUS OF TIRE MODELING - TIRE DESIGNER POINT OF VIEW 

The tire designer is interested in the load-carrying capability of the tire and 
hence requires a structural tire model for analytical studies. The development of 
reliable and cost effective structural tire models, however, is a formidable task. 
For example, the tire material properties are generally anisotropic and nonhomo- 
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geneous because of tire construction techniques. Such constituents as the nylon 
reinforcing cord and the rubber matrix tend to exhibit nonlinear elastic or visco- 
elastic behavior. Furthermore, typical loading conditions imposed upon the tire 
generally lead to large deformations and rotations and moderate strains. In 
addition, the material properties of the tire may be subject to temperature effects. 
The physics of tire loading and the thermal environment are not fully understood. 

All of these detrimental factors work in concert to complicate not only the modeling 
process but also the experimental techniques necessary to acquire corroborating data. 

Current tire modeling studies usually employ large general purpose finite 
element codes such as ADINA, MARC, STAGS, and NASTRAN as the modeling tool. 
Unfortunately these codes are not designed to handle the array of technical problems 
associated with tire modeling, and early tire modeling efforts tended to produce 
results which were below expectation and thus led to some skepticism within the 
U.S. tire industry. Moreover, there is no communication among the tire modelers 
within the tire industry because of proprietary concerns, and this lack of communi- 
cation can lead to a duplication of effort and the inability to apply a concerted 
effort to solve the problem. 

The consensus of the workshop participants was that the development of 
finite element technology for tire design is a nonlinear structural mechanics 
problem that should be addressed vigorously, and this effort should be 
directed by an agency outside the tire industry. It was recoimnended that this 
effort include three tasks: (1) establish a family of benchmark tire 

modeling problems, (2) produce a data base of experimental measurements x^hlch 
characterize the tire response to these problems, and (3) develop a number of 
special purpose computer codes that will model this response. 
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